*! Source of lmoremata11.mlib


*! {smcl}
*! {marker mm_density}{bf:mm_density.mata}{asis}
*! version 1.1.1  19nov2022  Ben Jann
version 11.2

// class & struct
local MAIN   mm_density
local SETUP  _`MAIN'_setup
local Setup  struct `SETUP' scalar
// real
local RS     real scalar
local RR     real rowvector
local RC     real colvector
local RV     real vector
// counters
local Int    real scalar
local IntC   real colvector
// string
local SS     string scalar
// boolean
local Bool   real scalar
local TRUE   1
local FALSE  0
// transmorphic
local T      transmorphic
local TS     transmorphic scalar
// pointers
local Pf     pointer(function) scalar
local PC     pointer(real colvector) scalar

mata:

// ---------------------------------------------------------------------------
// class definition
// ---------------------------------------------------------------------------

struct `SETUP' {
    // data
    `PC'    X             // pointer to X
    `PC'    w             // pointer to w
    `RS'    nobs          // number of obs/sum of weights
    `Bool'  pw            // weights are sampling weights
    `Bool'  sorted        // whether data is sorted

    // kernel
    `SS'    kernel        // name of kernel
    `Int'   adapt         // stages of adaptive estimator
    `Pf'    k             // kernel function
    `Pf'    kd            // kernel derivative
    `Pf'    kint          // kernel integral function
    `RS'    kh            // canonical bandwidth of kernel

    // bandwidth selection
    `RS'    h0            // user provided bandwidth
    `SS'    bwmethod      // bandwidth estimation method
    `RS'    adjust        // bandwidth adjustment factor
    `Int'   dpi           // number of DPI stages; default is 2
    `Bool'  qui           // quietly (omit SJPI/ISJ failure message)

    // support/boundary correction
    `RS'    lb            // lower boundary (missing if unbounded)
    `RS'    ub            // upper boundary (missing if unbounded)
    `SS'    bcmethod      // boundary correction method
    `Int'   bc            // 0=none, 1=renorm, 2=reflect, 3=linear correction
    `Bool'  rd            // relative data

    // other settings
    `Int'   n             // size of approximation grid
    `RS'    pad           // padding proportion of approximation grid
}

class `MAIN' {
    // settings
    public:
        void    data()        // set data
        `T'     kernel()      // kernel settings/retrieve kernel name
        `T'     bw()          // bandwidth settings
        `T'     support()     // support/boundary correction
        `T'     n()           // set/retrieve size of approximation grid
    private:
        void    new()         // initialize class with default settings
        void    clear()       // clear all results
        `Setup' setup         // settings
        `RC'    k(), kd(), kint() // raw kernel functions
        void    klc()         // compute terms for linear correction kernel
        `RC'    _K()          // kernel function including boundary correction
        void    checksuprt()  // check whether data is within support
    public:
        `RC'    X(), w()      // retrieve X and w
        `RS'    nobs()        // retrieve N (sum of weights)
        `Bool'  pw()          // retrieve pweighte flag
        `Bool'  sorted()      // retrieve sorted flag
        `Int'   adapt()       // retrieve stages of adaptive estimator
        `RS'    kh()          // retrieve canonical bandwidth of kernel
        `RC'    K()           // kernel function for external use
        `RC'    Kd()          // kernel derivative for external use
        `RS'    adjust()      // retrieve bw adjustment factor
        `Int'   dpi()         // retrieve dpi level
        `RS'    lb(), ub()    // retrieve lower and upper bounds of support
        `SS'    bc()          // retrieve boundary-correction method
        `Bool'  rd()          // retrieve relative data flag
        `RS'    pad()         // retrieve padding proportion
    
    // results
    public:
        `RC'    d()          // estimate/retrieve d
        `RS'    h()          // estimate/retrieve h
        `RC'    at()         // retrieve at
        `RC'    l()          // retrieve l
        `RC'    D()          // estimate/retrieve D
        `RC'    AT()         // set/retrieve AT
        `RC'    W()          // set/retrieve W
        `RC'    L()          // retrieve L
    private:
        `RS'    _h()         // estimate/retrieve h; return error if missing
        `RC'    d            // density estimate
        `RS'    h            // bandwidth
        `RC'    at           // evaluation grid
        `RC'    l            // local bandwidth factors (at observation level)
        `RC'    D            // full grid approximation estimate
        `RC'    D0           // adapt-1 approximation estimate
        `RC'    D0()         // adapt-1 approximation estimate
        `RC'    AT           // approximation grid
        `RC'    W            // grid counts
        `RC'    L            // local bandwidth factors
        
    // internal functions
    private:
        `RS'    h_sjpi(), h_isj(), h_dpi(), h_si(), h_ov(), h_no(), h_rd(), 
                _h_sjpi(), _h_isj(), h_root(), _h_root_fn()
        `Int'   _h_root()
        `RS'    df()
        `RC'    dd()
        void    dexact()
        `RC'    _dexact()
        `RC'    dapprox(), dapprox_fft(), dapprox_fft_rf(), dapprox_std()
        `RC'    lbwf()
        `RC'    ipolate()
        `RC'    grid()
        `RS'    scale()
}

// ---------------------------------------------------------------------------
// settings
// ---------------------------------------------------------------------------

void `MAIN'::new()
{
    kernel("")
    bw("")
    support(.)
    n(.)
}

void `MAIN'::clear()
{
    h  = .
    d = at = l = D = D0 = AT = W = L = J(0,1,.)
}

// D.data() -------------------------------------------------------------------

void `MAIN'::data(`RC' X, | `RC' w, `Bool' pw, `Bool' sorted)
{   // -sorted- indicates that
    //     X is sorted and non-missing
    //     w is non-missing and non-negative
    if (args()<2) w = 1
    if (args()<3) pw = `FALSE'
    if (args()<4) sorted = `FALSE'
    if (sorted==`FALSE') {
        if (missing(X) | missing(w)) _error(3351)
        if (any(w:<0)) {
            errprintf("{it:w} must not be negative\n")
            _error(3300)
        }
    }
    if (sorted==`FALSE') checksuprt(X, lb(), ub())
    setup.nobs   = mm_nobs(X, w)
    setup.X      = &X
    setup.w      = &w
    setup.pw     = (pw!=`FALSE')
    setup.sorted = (sorted!=`FALSE')
    clear()
}

`RC' `MAIN'::X()
{
    if (setup.X==NULL) return(J(0,1,.))
    return(*setup.X)
}

`RC' `MAIN'::w()
{
    if (setup.w==NULL) return(1)
    return(*setup.w)
}

`RS' `MAIN'::nobs() return(setup.nobs)

`Bool' `MAIN'::pw() return(setup.pw)

`Bool' `MAIN'::sorted() return(setup.sorted)

// D.kernel() -----------------------------------------------------------------

`T' `MAIN'::kernel(| `SS' kernel0, `Int' adapt)
{
    `SS' kernel
    
    // get
    if (args()==0) {
        return(setup.kernel)
    }
    // set
    if (adapt<0) _error(3300)
    kernel = strlower(strtrim(kernel0))
    if (kernel=="") kernel = "gaussian"  // default is "gaussian"
    setup.kernel = _mm_unabkern(strlower(strtrim(kernel)))
    setup.k      = _mm_findkern(setup.kernel)
    setup.kint   = _mm_findkint(setup.kernel)
    setup.kd     = _mm_findkderiv(setup.kernel)
    setup.kh     = (*_mm_findkdel0(setup.kernel))()
    setup.adapt = (adapt<. ? trunc(adapt) : 0) // default is 0
    clear()
}

`RS' `MAIN'::kh() return(setup.kh)

`Int' `MAIN'::adapt() return(setup.adapt)

`RC' `MAIN'::k(`RC' X) {
    return((*setup.k)(X))
}

`RC' `MAIN'::kd(`RC' X) {
    return((*setup.kd)(X))
}

`RC' `MAIN'::kint(`RS' l, | `RC' X)
{
    if (args()==1) return((*setup.kint)(l))
    return((*setup.kint)(l, X))
}

`RC' `MAIN'::K(`RC' X, `RC' x, `RC' h, | `Bool' fast)
{   // fast!=0: do not reset results to 0 outside support
    `RC' k
    
    k = _K(X, x, h):/h
    if (setup.bc==0)      return(k)
    if (fast & args()==4) return(k)
    if (lb()<.) k = k :* (X:>=lb() :& x:>=lb())
    if (ub()<.) k = k :* (X:<=ub() :& x:<=ub())
    return(k)
}

`RC' `MAIN'::_K(`RC' X, `RC' x, `RC' h)
{
    `RC' k, a0, a1, a2
    
    // unbounded support
    if (setup.bc==0) return(k((x:-X):/h))
    // renormalization
    if (setup.bc==1) {
        k = k((x:-X):/h)
        if      (ub()>=.) k = k :/  kint(1, (x:-lb()):/h) // lower bound only
        else if (lb()>=.) k = k :/  kint(1, (ub():-x):/h) // upper bound only
        else k = k :/ (kint(1, (ub():-x):/h) - kint(1, (lb():-x):/h))
        return(k)
    }
    // reflection
    if (setup.bc==2) {
        k = k((x:-X):/h)
        if (lb()<.) k = k + k((x:-2*lb():+X):/h)
        if (ub()<.) k = k + k((x:-2*ub():+X):/h)
        return(k)
    }
    // linear correction
    if (setup.bc==3) {
        klc(x, h, a0=., a1=., a2=.)
        k = (x:-X):/h
        return((a2 :- a1:*k):/(a0:*a2-a1:^2) :* k(k))
    }
    _error(3498) // not reached
}

void `MAIN'::klc(`RC' x, `RC' h, `RC' a0, `RC' a1, `RC' a2)
{
    `RC' l
    
    if (ub()>=.) { // lower bound only
        l  = (x:-lb()):/h
        a0 = kint(1, l)
        a1 = kint(3, l)
        a2 = kint(4, l)
        return
    }
    l = (ub():-x):/h
    a0 =  kint(1, l)
    a1 = -kint(3, l)
    a2 =  kint(4, l)
    if (lb()>=.) return // upper bound only
    l  = (lb():-x):/h
    a0 = a0 :- kint(1, l)
    a1 = a1 :+ kint(3, l)
    a2 = a2 :- kint(4, l)
}

`RC' `MAIN'::Kd(`RC' X, `RC' x, `RC' h, | `Bool' fast)
{   // fast!=0: do not reset results to 0 outside support
    `RC' k, a0, a1, a2
    
    // unbounded support
    if (setup.bc==0) return(-kd((x:-X):/h) :/ h:^2)
    // renormalization
    if (setup.bc==1) {
        k = kd((x:-X):/h)
        if      (ub()>=.) k = k :/  kint(1, (x:-lb()):/h) // lower bound only
        else if (lb()>=.) k = k :/  kint(1, (ub():-x):/h) // upper bound only
        else k = k :/ (kint(1, (ub():-x):/h) - kint(1, (lb():-x):/h))
        k = -k :/ h:^2
    }
    // reflection
    else if (setup.bc==2) {
        k = kd((x:-X):/h)
        if (lb()<.) k = k - kd((x:-2*lb():+X):/h)
        if (ub()<.) k = k - kd((x:-2*ub():+X):/h)
        k = -k :/ h:^2
    }
    // linear correction
    else if (setup.bc==3) {
        klc(x, h, a0=., a1=., a2=.)
        k = (x:-X):/h
        k = (a1:*k(k) - (a2 :- a1:*k):*kd(k)) :/ ((a0:*a2-a1:^2):*h:^2)
    }
    else _error(3498) // not reached
    // reset to zero outside of support
    if (fast & args()==4) return(k)
    if (lb()<.) k = k :* (X:>=lb() :& x:>=lb())
    if (ub()<.) k = k :* (X:<=ub() :& x:<=ub())
    return(k)
}

// D.bw() -------------------------------------------------------------------

`T' `MAIN'::bw(| `TS' bw, `RS' adj, `Int' dpi, `Bool' qui)
{
    `RS' h
    `SS' method
    
    // get
    if (args()==0) {
        if (setup.h0<.) return(setup.h0)
        return(setup.bwmethod)
    }
    // set
    if (adj<=0) _error(3300)
    if (dpi<0)  _error(3300)
    if (args()<4) qui = `FALSE'
    if (!isstring(bw)) {
        if (bw<=0) _error(3300)
        if (bw<.) h = bw
    }
    else method = bw
    method = mm_strexpand(strlower(strtrim(method)),
        ("silverman", "normalscale", "oversmoothed", "sjpi", "dpi", "isj"),
        "sjpi") // default is "sjpi"
    setup.h0       = h
    setup.bwmethod = method
    setup.adjust   = (adj<. ? adj : 1)          // default is 1
    setup.dpi      = (dpi<. ? trunc(dpi) : 2)   // default is 2
    setup.qui      = (qui!=`FALSE')
    clear()
}

`RS' `MAIN'::adjust() return(setup.adjust)

`Int' `MAIN'::dpi() return(setup.dpi)

// D.support() ----------------------------------------------------------------

`T' `MAIN'::support(| `RV' minmax, `SS' method, `Bool' rd)
{
    `RS' lb, ub
    
    // get
    if (args()==0) {
        return((setup.lb, setup.ub))
    }
    // set
    if (args()<3) rd = `FALSE'
    if (length(minmax)>2) _error(3200)
    if (length(minmax))   lb = minmax[1]
    else                  lb = .
    if (length(minmax)>1) ub = minmax[2]
    else                  ub = .
    if (lb<. & ub<.) {
        if (lb>=ub) _error(3300)
    }
    if (rd) {
        if (lb >= .)  lb = 0
        if (lb <  0)  _error(3300)
        if (ub >= .)  ub = 1
        if (ub >  1)  _error(3300)
        if (lb > ub)  _error(3300)
    }
    if (sorted()==`FALSE') checksuprt(X(), lb, ub)
    setup.bcmethod = mm_strexpand(strlower(strtrim(stritrim(method))),
        ("renormalization", "reflection", "linear correction"),
        "renormalization") // default is "renormalization"
    setup.lb = lb
    setup.ub = ub
    setup.rd = (rd!=`FALSE')
    if (setup.lb>=. & setup.ub>=.)                  setup.bc = 0
    else if (setup.bcmethod=="renormalization")     setup.bc = 1
    else if (setup.bcmethod=="reflection")          setup.bc = 2
    else if (setup.bcmethod=="linear correction")   setup.bc = 3
    else _error(3498) // cannot be reached
    clear()
}

`RS' `MAIN'::lb() return(setup.lb)

`RS' `MAIN'::ub() return(setup.ub)

`SS' `MAIN'::bc() return(setup.bcmethod)

`Bool' `MAIN'::rd() return(setup.rd)

void `MAIN'::checksuprt(`RC' X, `RS' lb, `RS' ub)
{
    if (rows(X)==0) return // no data
    if (lb<.) {
        if (min(X)<lb) {
            errprintf("{it:X} contains values out of support\n")
            _error(3300)
        }
    }
    if (ub<.) {
        if (max(X)>ub) {
            errprintf("{it:X} contains values out of support\n")
            _error(3300)
        }
    }
}

// D.n() ----------------------------------------------------------------------

`T' `MAIN'::n(| `Int' n, `RS' pad)
{
    // get
    if (args()==0) return(setup.n)
    // set
    if (n<1)   _error(3300)
    if (pad<0) _error(3300)
    setup.n   = (n<. ? trunc(n) : 1024)  // default is 1024
    setup.pad = (pad<. ? pad : 0.1)      // default is 0.1
    clear()
}

`RS' `MAIN'::pad() return(setup.pad)

// ---------------------------------------------------------------------------
// bandwidth selection
// ---------------------------------------------------------------------------

`RS' `MAIN'::h() {
    if (h<.) return(h)
    // user bandwidth
    if (setup.h0<.) { 
        h = setup.h0 * adjust()
        return(h)
    }
    // data-driven bandwidth selection
    if (bw()=="sjpi") {
        h = h_sjpi()
        if (h>=.) {
            if (setup.qui==`FALSE')
                display("{txt}(SJPI bandwidth estimation failed; using DPI method)")
            h = h_dpi()
        }
    }
    else if (bw()=="isj") {
        h = h_isj()
        if (h>=.) {
            if (setup.qui==`FALSE')
                display("{txt}(ISJ bandwidth estimation failed; using DPI method)")
            h = h_dpi()
        }
    }
    else if (bw()=="dpi")          h = h_dpi()
    else if (bw()=="silverman")    h = h_si()
    else if (bw()=="oversmoothed") h = h_ov()
    else if (bw()=="normalscale")  h = h_no()
    else _error(3498) // cannot be reached
    // correction for pweights
    if (pw() & rows(w())!=1) {
        h = h * (sum(w():^2)/(nobs()))^.2
    }
    // final adjustments
    h = h * kh() * adjust()
    if (h<=0) h = . // h must be strictly positive
    return(h)
}

`RS' `MAIN'::h_rd(`RS' s) // relative data correction factor
{
    if (rd()) return((1 + 1/(2 * sqrt(pi()) * s))^.2)
    return(1)
}

`RS' `MAIN'::_h()
{
    if (h()>=.) {
        errprintf("bandwidth could not be determined\n")
        _error(3498)
    }
    return(h)
}

// Sheather-Jones solve-the-equation plug-in selection rule -------------------

`RS' `MAIN'::h_sjpi()
{
    `RS'  n, s, hmin, h_os, sda, tdb, tdc, alpha, beta
    `RC'  AT, W
    
    AT = AT()
    W = W()
    s = scale(0, AT, W, 1)           // min of sd and iqr
    if (pw()) {                      // pweights: normalize grid counts
        n = rows(X())
        W = W * (n / nobs())
    }
    else n = nobs()
    sda = df(AT, W, n, 1.241 * s * n^(-1/7), 2)
    tdb = df(AT, W, n, 1.230 * s * n^(-1/9), 3)
    alpha = 1.357 * (sda/tdb)^(1/7)
    if (rd()) {
        tdc = df(AT, W, n, 1.304 * s * n^(-1/5), 1)
        beta = 1.414 * (sda/tdc)^(1/3)
    }
    hmin = .5 * (AT[n()]-AT[1])/(n()-1) * mm_kdel0_gaussian()/mm_kdel0_rectangle()
    h_os = s * (243/(35*n))^.2 * mm_kdel0_gaussian() * h_rd(s) // oversmoothed (modified)
        // // plot objective function:
        // `Int' i
        // `RC'  at, hh
        // at = rangen(hmin, 3*h_os, 100)
        // hh = J(rows(at),1,.)
        // for (i=rows(at);i;i--) hh[i] = _h_sjpi(at[i], AT, W, n, alpha, beta)
        // mm_plot((hh,at),"line",sprintf("xline(%g)",h_os))
    return(h_root(0, hmin, h_os, AT, W, n, alpha, beta) * 
        (n/nobs())^.2 / mm_kdel0_gaussian())
}

`RS' `MAIN'::h_root(`Int' m, `RS' hmin, `RS' h_os,
    | `T' o1, `T' o2, `T' o3, `T' o4, `T' o5)
{
    `RS'  ax, bx, h
    `Int' rc
    
    bx = h_os
    ax = max((hmin, bx/2))
    rc = _h_root(m, h=., ax, bx, o1, o2, o3, o4, o5)
    if (rc==2) { // move down
        bx = ax
        while (1) {
            if (bx<=hmin) break // cannot go below hmin
            ax = max((hmin, bx/2))
            rc = _h_root(m, h, ax, bx, o1, o2, o3, o4, o5)
            if (rc==2) bx = ax // continue moving down
            else break // this also stops if there is a change in 
                      //  direction (rc==3) (in this case: h = current bx)
        }
    }
    else if (rc==3) { // move up
        ax = bx; bx = ax*1.5
        while (1) {
            rc = _h_root(m, h, ax, bx, o1, o2, o3, o4, o5)
            if (rc==3) {; ax = bx; bx = ax*1.5; } // continue moving up
            else break // this also stops if there is a change in 
                       //  direction (rc==2) (in this case: h = current ax)
        }
    }
    if (h<=hmin) { // solution is smaller than hmin
        ax = h_os; bx = ax*1.5
        rc = _h_root(m, h, ax, bx, o1, o2, o3, o4, o5)
        if (rc==2) return(.) // moving up does not help
        if (rc==3) { // continue moving up
            ax = bx; bx = ax*1.5
            while (1) {
                rc = _h_root(m, h, ax, bx, o1, o2, o3, o4, o5)
                if (rc==3) {; ax = bx; bx = ax*1.5; } // continue moving up
                else break // this also stops if there is a change in 
                           //  direction (rc==2) (in this case: h = current ax)
            }
        }
    }
    return(h)
}

`Int' `MAIN'::_h_root(`Int' m, `RS' x, `RS' ax, `RS' bx, 
    | `T' o1, `T' o2, `T' o3, `T' o4, `T' o5)
{   // root finder for SJPI (m!=1) and ISJ (m=1); adapted from mm_root()
    `Int' maxit, itr
    `RS'  tol, tol_act
    `RS'  a, b, c, fa, fb, fc, prev_step, p, q, new_step, t1, cb, t2

    tol   = 0
    maxit = 100
    x = .; a = ax; b = bx
    fa = _h_root_fn(m, a, o1, o2, o3, o4, o5)
    fb = _h_root_fn(m, b, o1, o2, o3, o4, o5)
    c = a; fc = fa
    if ( fa==. ) return(0) // abort if fa missing => x=.
    if ( (fa > 0 & fb > 0) | (fa < 0 & fb < 0) ) {
        if ( abs(fa) < abs(fb) ) {
            x = a; return(2)
        }
        x = b; return(3)
    }
    for (itr=1; itr<=maxit; itr++) {
        if ( fb==. ) return(0)
        prev_step = b-a
        if ( abs(fc) < abs(fb) ) {
            a = b;  b = c;  c = a; fa = fb; fb = fc; fc = fa
        }
        tol_act = 2*epsilon(b) + tol/2
        new_step = (c-b)/2
        if ( abs(new_step) <= tol_act | fb == 0 ) {
             x = b
             return(0)
        }
        if ( abs(prev_step) >= tol_act & abs(fa) > abs(fb) ) {
            cb = c-b
            if ( a==c ) {
                t1 = fb/fa
                p = cb*t1
                q = 1.0 - t1
            }
            else {
                q = fa/fc;  t1 = fb/fc;  t2 = fb/fa
                p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) )
                q = (q-1.0) * (t1-1.0) * (t2-1.0)
            }
            if ( p>0 ) q = -q
            else      p = -p
            if ( p < (0.75*cb*q-abs(tol_act*q)/2) & p < abs(prev_step*q/2) )
                new_step = p/q
        }
        if ( abs(new_step) < tol_act ) {
            if ( new_step > 0 ) new_step = tol_act
            else                new_step = -tol_act
        }
        a = b;  fa = fb
        b = b + new_step; fb = _h_root_fn(m, b, o1, o2, o3, o4, o5)
        if ( (fb > 0 & fc > 0) | (fb < 0 & fc < 0) ) {
            c = a;  fc = fa
        }
    }
    x = b
    return(1)
}

`RS' `MAIN'::_h_root_fn(`Int' m, `RS' x, | `T' o1, `T' o2, `T' o3, `T' o4, `T' o5)
{
    if (m==1) return(_h_isj(x, o1, o2, o3))   // ISJ method
    return(_h_sjpi(x, o1, o2, o3, o4, o5))    // SJPI method
}

`RS' `MAIN'::_h_sjpi(`RS' h, `RC' AT, `RC' W, `RS' n, `RS' alpha, `RS' beta)
{
    `RS' d
    
    d = mm_kint_gaussian(2)
    if (rd()) d = d * (1 + df(AT, W, n, beta * h^(5/3), 0))
    return((d / (n * df(AT, W, n, alpha * h^(5/7), 2)))^0.2 - h)
}

`RS' `MAIN'::df(`RC' AT, `RC' W, `RS' n, `RS' h, `Int' d)
{
    return( (-1)^d * sum(W :* dd(AT, W, n, h, d)) / n )
}

`RC' `MAIN'::dd(`RC' AT, `RC' W, `RS' n, `RS' h, `Int' d)
{    // d'th density derivative using gaussian kernel
    `Int' M, L, i, first, last
    `RS'  a, b
    `RC'  kappam, arg, hmold0, hmold1, hmnew, w
    
    // compute kappam
    M = rows(AT)
    a = AT[1]; b = AT[M]
    L = (M-1) * (1 + lb()<. + ub()<.)
    L = max( (min( (floor((4 + 2*d) * h * (M-1)/(b-a)), L) ), 1) )
    arg = (0::L) * (b-a) / (h*(M-1))
    kappam = normalden(arg)
    hmold0 = 1; hmold1 = arg; hmnew  = 1
    for (i=2; i<=2*d; i++) { // compute mth degree Hermite polynomial
        hmnew  = arg :* hmold1 :- (i-1) * hmold0
        hmold0 = hmold1; hmold1 = hmnew
    }
    kappam = hmnew :* kappam
    
    // unbounded estimator
    if (lb()>=. & lb()>=.) {
        return( convolve((kappam[L+1::1] \ kappam[|2 \ L+1|]), 
            W)[|L+1 \ L+M|] / (n*h^(2*d+1)) )
    }
    
    // reflection estimator
    w = W
    if (lb()<.) {
        first = M
        w = W[M::2] \ w
        w[first] = 2 * w[first]
    }
    else first = 1
    last = first + (M-1)
    if (ub()<.) {
        w = w \ W[M-1::1]
        w[last] = 2 * w[last]
    }
    return( convolve((kappam[L+1::1] \ kappam[|2 \ L+1|]), 
        w)[|L+first \ L+last|] / (n*h^(2*d+1)) )
}

// "improved" SJPI (diffusion method) -----------------------------------------

`RS' `MAIN'::h_isj()
{
    `Int' n
    `RS'  N, s, hmin, h_os
    `RC'  AT, W, a
    
    // step 1: bin data on regular grid
    n = 2^ceil(ln(n())/ln(2))        // round up to next power of 2
    n = max((n, 1024))               // enforce min grid size of at least 1024
    AT = grid(n)                     // generate grid
    if (sorted()) W = _mm_exactbin(X(), w(), grid(n+1))
    else          W = mm_fastexactbin(X(), w(), grid(n+1))
        // need to use exact binning because linear binning would introduce
        // some (non-vanishing) bias at the boundaries (doubling the first and
        // last grid count does not seem to help); a consequence of exact 
        // binning is that the density estimate will be slightly shifted/stretched
        // to the left; this error can be substantial if the grid size is small,
        // but it vanishes with increasing grid size
    s = scale(0, AT, W, 1)           // min of sd and iqr
    W = W / nobs()                   // relative frequencies
    if (pw()) N = rows(X())          // obtain sample size
    else      N = nobs()
    // step 2: obtain discrete cosine transform of binned data
    a = Re( (1 \ 2 * exp(1i * (1::n-1) * pi() / (2*n)))
         :* fft(W[mm_seq(1,n-1,2)] \ W[mm_seq(n,2,2)]) )
    // step 3: compute bandwidth
    hmin = (.5/(n-1) * mm_kdel0_gaussian()/mm_kdel0_rectangle())^2
    h_os = (s * (243/(35*N))^.2 * mm_kdel0_gaussian() * h_rd(s) / (AT[n]-AT[1]))^2
        // // plot objective function:
        // `Int' i
        // `RC'  at, hh
        // at = rangen(hmin, 3*h_os, 100)
        // hh = J(rows(at),1,.)
        // for (i=rows(at);i;i--) hh[i] = _h_isj(at[i], N, (1::n-1):^2, (a[2::n]/2):^2)
        // mm_plot((hh,at),"line",sprintf("xline(%g)",h_os))
    return(sqrt(h_root(1, hmin, h_os, N, (1::n-1):^2, (a[2::n]/2):^2)) * 
        (AT[n]-AT[1]) * (N/nobs())^.2 / mm_kdel0_gaussian() * h_rd(s))
}

`RS' `MAIN'::_h_isj(`RS' h, `RS' N, `RC' I, `RC' a2)
{
    `Int' l, s
    `RS'  K0, c
    `RC'  f, t
    
    l = 7
    f = 2 * pi()^(2*l) * sum(I:^l :* a2 :* exp(-I * pi()^2 * h))
    for (s=l-1; s>=2; s--) {
        K0 = mm_prod(mm_seq(1, 2*s-1, 2)) / sqrt(2*pi())
        c  = (1 + (1/2)^(s + 1/2)) / 3
        t  = (2 * c * K0/N/f):^(2/(3 + 2*s))
        f  = 2 * pi()^(2*s) * sum(I:^s :* a2 :* exp(-I * pi()^2 * t))
    }
    return((2 * N * sqrt(pi()) * f)^(-2/5) - h)
}

// Sheather-Jones direct plug-in selection rule -------------------------------

`RS' `MAIN'::h_dpi()
{
    `RS'  n, s, alpha, psi, psi0, alpha0
    `Int' i
    `RC'  AT, W
    
    i = dpi()
    if (i==0) return(h_no()) // h normalscale
    else {
        AT = AT()
        W  = W()
        s = scale(0, AT, W, 1)           // min of sd and iqr
        if (pw()) {                      // pweights: normalize grid counts
            n = rows(X())
            W = W * (n / nobs())
        }
        else n = nobs()
        alpha = (2 * (sqrt(2) * s)^(3 + 2 * (i+1)) /
                ((1 + 2 * (i+1)) * n))^(1/(3 + 2 * (i+1)))
        if (rd()) alpha0 = (2 * (sqrt(2) * s)^(3 + 2 * (i-1)) /
                           ((1 + 2 * (i-1)) * n))^(1/(3 + 2 * (i-1)))
        for (; i; i--) {
            psi = df(AT, W, n, alpha, i+1)
            if (rd()) psi0 = df(AT, W, n, alpha0, i-1)
            else      psi0 = 0
            if (i>1) {
                alpha = ( factorial(i*2) / (2^i * factorial(i)) *
                          sqrt(2/pi()) / (psi*n) )^(1/(3 + 2*i))
                if (rd()) alpha0 = ( factorial((i-2)*2) / (2^(i-2) * 
                                   factorial(i-2)) * sqrt(2/pi()) / 
                                   (psi0*n) )^(1/(3 + 2*(i-2)))
            }
        }
    }
    return( ((1 + psi0) / (psi * nobs()))^.2 )
}

// optimal of Silverman selection rule ----------------------------------------

`RS' `MAIN'::h_si()
{
    `RS' s
    
    s = scale(0, X(), w(), sorted()) // min of sd and iqr
    return(0.9/mm_kdel0_gaussian() * s / nobs()^.2 * h_rd(s))
}

// oversmoothed selection rule ------------------------------------------------

`RS' `MAIN'::h_ov()
{
    `RS' s
    
    s = scale(1, X(), w(), sorted()) // sd
    return((243/35)^.2 * s / nobs()^.2 * h_rd(s))
}

// normal scale selection rule ------------------------------------------------

`RS' `MAIN'::h_no()
{
    `RS' s
    
    s = scale(0, X(), w(), sorted()) // min of sd and iqr
    return((8*sqrt(pi())/3)^.2 * s / nobs()^.2 * h_rd(s))
}

// ---------------------------------------------------------------------------
// density estimation
// ---------------------------------------------------------------------------

`RC' `MAIN'::d(| `RV' o1, `RS' o2, `RS' o3, `RS' o4)
{
    // case 0: return existing result
    if (args()==0) return(d)
    // case 1: o1 contains grid
    if (args()<=2) {
        if (cols(o1)!=1)   at = o1'
        else               at = o1
        // case 1a: exact estimator
        if (args()==2 & o2) dexact()
        // case 1b: approximation estimator
        else d = ipolate(AT(), D(), at, mm_issorted(at))
        return(d)
    }
    // case 2: o1=n, o2=from, o3=to
    at = grid(o1, _h(), o2, o3)
    // case 2a: exact estimator
    if (args()==4 & o4) dexact()
    // case 2b: approximation estimator
    else d = ipolate(AT(), D(), at, mm_issorted(at))
    return(d)
}

`RC' `MAIN'::at() return(at)

// exact estimator ------------------------------------------------------------

void `MAIN'::dexact()
{
    `Int'  n
    `IntC' p
    
    n = rows(at)
    if (setup.bc & n>0) {
        // check for evaluation points out of support and set density
        // to zero for these points
        if (lb()<. & ub()<.) p = select(1::n, at:>=lb() :& at:<=ub())
        else if (lb()<.)     p = select(1::n, at:>=lb())
        else if (ub()<.)     p = select(1::n, at:<=ub())
        if (length(p)!=n) {
            d = J(n,1,0)
            if (length(p)) d[p] = _dexact(X(), w(), _h() :* l(), at[p])
            return
        }
    }
    d = _dexact(X(), w(), _h() :* l(), at)
}

`RC' `MAIN'::_dexact(`RC' x, `RC' w, `RC' h, `RC' at)
{
    `Int' i
    `RC'  d
    
    i = rows(at)
    d = J(i,1,.)
    // using slightly different method depending on whether h and w are
    // scalar or not to save a bit of computer time, if possible
    if (rows(h)==1 & rows(w)==1) {
        for (;i;i--) d[i] = w/h * sum(_K(x, at[i], h))
    }
    else if (rows(h)==1) {
        for (;i;i--) d[i] = sum(w :* _K(x, at[i], h)) / h
    }
    else if (rows(w)==1) {
        for (;i;i--) d[i] = w * sum(_K(x, at[i], h) :/ h)
    }
    else {
        for (;i;i--) d[i] = sum(w:/h :* _K(x, at[i], h))
    }
    return(d / nobs())
}

`RC' `MAIN'::l()
{
    if (rows(l)) return(l)
    if (!adapt()) return(1)
    return(lbwf(ipolate(AT(), D0(), X(), sorted()), w()))
}

// binned approximation estimator ---------------------------------------------

`RC' `MAIN'::D()
{
    if (rows(D)) return(D)
    if (adapt()) (void) D0()
    D = dapprox()
    return(D)
}

`RC' `MAIN'::D0()
{
    `Int' i
    
    if (rows(D0)) return(D0)
    for (i=adapt();i;i--) D0 = dapprox()
    return(D0)
}

`RC' `MAIN'::AT()
{
    if (rows(AT)) return(AT)
    AT = grid(n())
    return(AT)
}

`RC' `MAIN'::W()
{
    if (rows(W)) return(W)
    if (sorted()) W = _mm_linbin(X(), w(), AT())
    else          W = mm_fastlinbin(X(), w(), AT())
    return(W)
}

`RC' `MAIN'::L()
{
    if (rows(L)) return(L)
    if (!adapt()) return(1)
    L = lbwf(D0(), W())
    return(L)
}

`RC' `MAIN'::dapprox()
{
    `RC' h
    
    // create grid and compute grid counts if necessary
    if (rows(W)==0) (void) W()
    // obtain h
    h = _h()
    if (rows(D0)) h = h :* lbwf(D0, W) // adaptive estimator
    // FFT estimation if h is constant
    if (rows(h)==1) {
        if (setup.bc==3) return(dapprox_std(h))    // linear correction (no FFT)
        if (setup.bc==2) return(dapprox_fft_rf(h)) // reflection
                         return(dapprox_fft(h))    // no bc or renormalization
    }
    // else use standard estimator
    return(dapprox_std(h))
}

`RC' `MAIN'::dapprox_fft(`RS' h)
{
    `Int' M, L
    `RS'  a, b, tau
    `RC'  kappa, D
    
    M = rows(AT)
    a = AT[1]; b = AT[M]
    L = M - 1
    // reduce number of evaluation points if possible
    if (kernel()!="gaussian") {
        if (kernel()=="cosine")            tau = .5
        else if (kernel()=="epanechnikov") tau = sqrt(5)
        else                               tau = 1
        L = max( (min( (floor(tau*h*(M-1)/(b-a)), L) ), 1) )
    }
    // compute kappa and obtain FFT
    kappa = k( (0::L) * (b-a) / (h*(M-1)) )
    D = convolve((kappa[L+1::1]\kappa[|2 \ L+1|]), W)[|L+1 \ L+M|] / (nobs()*h)
    if (setup.bc==0) return(D)
    // renormalization boundary correction
    if (ub()>=.) return(D :/  kint(1, (AT:-lb()):/h)) // lower bound only
    if (lb()>=.) return(D :/  kint(1, (ub():-AT):/h)) // upper bound only
    return(D :/ (kint(1, (ub():-AT):/h) - kint(1, (lb():-AT):/h)))
}

`RC' `MAIN'::dapprox_fft_rf(`RS' h)
{
    `Int' M, L, first, last
    `RS'  a, b, tau
    `RC'  kappa, w
    
    M = rows(AT)
    a = AT[1]; b = AT[M]
    L = (M-1) * (1 + (lb()<.) + (ub()<.))
    // reduce number of evaluation points if possible
    if (kernel()!="gaussian") {
        if (kernel()=="cosine")            tau = .5
        else if (kernel()=="epanechnikov") tau = sqrt(5)
        else                               tau = 1
        L = max( (min( (floor(tau*h*(M-1)/(b-a)), L) ), 1) )
    }
    // expand vector of grid counts
    w = W
    if (lb()<.) {
        first = M
        w = W[M::2] \ w
        w[first] = 2 * w[first]
    }
    else first = 1
    last = first + (M-1)
    if (ub()<.) {
        w = w \ W[M-1::1]
        w[last] = 2 * w[last]
    }
    // compute kappa and obtain FFT
    kappa = k( (0::L) * (b-a) / (h*(M-1)) )
    return(convolve((kappa[L+1::1]\kappa[|2 \ L+1|]), w)[|L+first \ L+last|] / 
        (nobs()*h))
}

`RC' `MAIN'::dapprox_std(`RC' h)
{
    `Int' n, i, a, b
    `RC'  r, D
    
    // no computational shortcut in case of gaussian kernel
    if (kernel()=="gaussian") return(_dexact(AT, W, h, AT))
    // other kernels: restrict computation to relevant range of evaluation points
    n = n()
    r = h
    if (kernel()=="cosine")            r = r * .5
    else if (kernel()=="epanechnikov") r = r * sqrt(5)
    r = trunc(r * (n-1) / (AT[n]-AT[1])) :+ 1 // add 1 to prevent roundoff error
    D = J(n,1,0)
    if (rows(h)==1) {
        for (i=n;i;i--) {
            a = max((1, i-r))
            b = min((n, i+r))
            D[|a \ b|] = D[|a \ b|] + W[i] / h * _K(AT[i], AT[|a \ b|], h)
        }
        return(D :/ nobs())
    }
    for (i=n;i;i--) {
        a = max((1, i-r[i]))
        b = min((n, i+r[i]))
        D[|a \ b|] = D[|a \ b|] + W[i] / h[i] * _K(AT[i], AT[|a \ b|], h[i])
    }
    return(D :/ nobs())
}

// ---------------------------------------------------------------------------
// helper functions
// ---------------------------------------------------------------------------

`RC' `MAIN'::grid(`Int' n, | `RS' h, `RS' from, `RS' to)
{
    `RS' tau
    `RR' range, minmax
    
    range = J(1,2,.)
    if (from<.)                 range[1] = from
    else if (lb()<. & from<.y)  range[1] = lb()
    if (to<.)                   range[2] = to
    else if (ub()<. & to<.y)    range[2] = ub()
    if (missing(range)) {
        minmax = minmax(X())
        if (range[1]>=. & from==.z) range[1] = minmax[1]
        if (range[2]>=. & to==.z)   range[2] = minmax[2]
        if (missing(range)) {
            // if h is not provided:
            // - extend grid below min(x) and above max(x) by pad()% of data range
            // if h is provided:
            // - extend grid by the kernel halfwidth such that density can go to 
            //   zero outside data range (only approximately for gaussian kernel or
            //   in case of adaptive estimator), but limit by pad()% of data range
            tau = (minmax[2]-minmax[1]) * pad()
            if (h<.) tau = min((tau, h * (kernel()=="epanechnikov" ? sqrt(5) : 
                (kernel()=="cosine" ? .5 : (kernel()=="gaussian" ? 3 : 1)))))
            minmax[1] = minmax[1] - tau
            minmax[2] = minmax[2] + tau
            if (range[1]>=.) {
                if (lb()<.) minmax[1] = max((lb(),minmax[1]))
                range[1] = minmax[1]
            }
            if (range[2]>=.) {
                if (ub()<.) minmax[2] = min((ub(),minmax[2]))
                range[2] = minmax[2]
            }
        }
    }
    if (range[1]>range[2]) _error(3300)
    return(rangen(range[1], range[2], n))
}

`RS' `MAIN'::scale(`Int' type, `RC' X, `RC' w, `Bool' sorted)
{   // type: 0 = min(sd,iqr), 1 = sd, 2 = iqr
    // iqr will be replaced by sd if 0
    `RS' iqr, sd
    
    if (type!=1) {
        if (sorted) iqr = _mm_iqrange(X, w) / 1.349
        else        iqr =  mm_iqrange(X, w) / 1.349
    }
    if (type!=2 | iqr<=0) {
        sd = sqrt(variance(X, w))
        if (pw()) sd = sd * sqrt( (nobs()-1) / (nobs() - nobs()/rows(X())) )
    }
    if (type==1) return(sd)
    if (iqr<=0)  iqr = sd
    if (type==2) return(iqr)
    return(min((sd, iqr)))
}

`RC' `MAIN'::lbwf(`RC' d, `RC' w) // local bandwidth factors
{
    `RC' l
    
    l = sqrt( exp(mean(log(d), w)) :/ d)  // exp(...) -> geometric mean
    return(editmissing(l, 1))
}

`RC' `MAIN'::ipolate(`RC' AT, `RC' D, `RC' at, `Bool' sorted)
{
    `RC' d, p
    
    if (sorted) d = mm_fastipolate(AT, D, at)
    else {
        p = order(at, 1)
        d = mm_fastipolate(AT, D, at[p])
        d[p] = d
    }
    _editmissing(d, 0)  // set density outside of grid to 0
    return(d)
}

end



*! {smcl}
*! {marker mm_qr}{bf:mm_qr.mata}{asis}
*! version 1.0.5  30mar2021  Ben Jann
version 11
mata:

class mm_qr
{
    // constructors
    private:
        void           new()
        void           init() // default settings
        void           clear(), clear1(), clear2(), clear2b(), clear3()

    // setup
    public:
        void           data()
        transmorphic   qd()
        transmorphic   demean()
        transmorphic   collin()
        transmorphic   p()
        transmorphic   b_init()
        transmorphic   tol()
        transmorphic   maxiter()
        transmorphic   beta()
        transmorphic   method()
        transmorphic   log()
    
    private:
        pointer scalar y
        pointer scalar X
        pointer scalar w
        real scalar    n, N
        real scalar    cons
        real scalar    k, kadj
        real scalar    K, Kadj
        real scalar    p
        real colvector b_init, b_ls
        real scalar    b_init_user
        real scalar    maxiter
        real scalar    tol
        real scalar    beta
        string scalar  method
        real scalar    qd
        real scalar    demean
        real scalar    collin
        real colvector omit
        real rowvector xindx
        real scalar    log
    
    // results
    public:
        real colvector b()
        real colvector xb()
        real scalar    gap()
        real scalar    sdev()
        real scalar    iter()
        real scalar    converged()
        real scalar    n() // number of observations
        real scalar    N() // sum of weights
        real scalar    cons()
        real scalar    k() // number of predictors
        real scalar    K() // number of coefficients
        real colvector omit()
        real scalar    k_omit()
        real scalar    ymean()
        real rowvector means()
    
    private:
        real colvector b
        real scalar    conv
        real scalar    iter
        real scalar    gap
        real scalar    sdev
        real scalar    ymean
        real rowvector means
        
    // functions
    private:
        void           err_nodata()
        void           printflush()
        void           printlog()
        real colvector _xb()
        real scalar    _sdev()
        void           set_collin()
        transmorphic   lsfit()
        void           _b_init()
        real colvector rmomit(), addomit()
        real colvector meanadj()
        void           fit(), fnb()
        pointer scalar get_X()
        void           gen_z_w()
        real matrix    cross()
        real scalar    minselect()
}

void mm_qr::new()
{
    init()
    clear()
}

void mm_qr::init()
{
    p         = .5
    qd        = 1
    demean    = 1
    collin    = 1
    tol       = 1e-8
    maxiter   = st_numscalar("c(maxiter)")
    beta      = 0.99995
    method    = "fnb"
    log       = 0
}

void mm_qr::clear() // if new() or data()
{
    y = X = w = NULL
    n = N = cons = k = K = ymean = means = .z
    b_init_user = 0 // b_init can only be set after data has been set
    clear1()
}

void mm_qr::clear1() // if qd() or demean()
{
    b_ls = .z
    if (!b_init_user) b_init = .z
    clear2()
}

void mm_qr::clear2() // if collin()
{
    kadj = Kadj = omit = xindx = .z
    clear3()
}

void mm_qr::clear2b() // if p()
{
    if (!b_init_user & cons) b_init = .z
    clear3()
}

void mm_qr::clear3() // if b_init(), tol(), maxiter(), beta()
{
    b = conv = iter = gap = sdev = .z
}

void mm_qr::err_nodata()
{
    if (y==NULL) {
        display("{err}data not set")
        _error(3498)
    }
}

void mm_qr::printflush(string scalar s)
{
    printf(s)
    displayflush()
}

void mm_qr::printlog(real scalar d, real scalar iter)
{
    printf("{txt}")
    printf("{txt}Iteration %g:", iter)
    printf("{col 16}mreldif() in b = {res}%11.0g;", d)
    printf("{txt}  duality gap = {res}%11.0g\n", gap)
    displayflush()
}

void mm_qr::data(real colvector y0, | real matrix X0, real colvector w0,
    real colvector cons0)
{
    clear()
    // do error checks before storing anything
    if (!(X0==. | X0==J(0,0,.))) {
        if (rows(X0)!=rows(y0)) _error(3200)
    }
    if (rows(w0)) {
        if (rows(w0)!=rows(y0) & rows(w0)!=1) _error(3200)
    }
    // now start storing
    y = &y0
    n = rows(y0)
    if (X0==. | X0==J(0,0,.)) X = &J(n, 0, 1)
    else                      X = &X0
    if (rows(w0)) {
        if (rows(w0)==1) N = n * w0
        else             N = quadsum(w0)
        w = &w0
    }
    else {
        N = n
        w = &1
    }
    cons = (cons0!=0)
    k = cols(*X)
    K = k + cons
}

transmorphic mm_qr::qd(| real scalar qd0)
{
    if (args()==0) return(qd)
    if (qd==(qd0!=0)) return // no change
    qd = (qd0!=0)
    clear1()
}

transmorphic mm_qr::demean(| real scalar demean0)
{
    if (args()==0) return(demean)
    if (demean==(demean0!=0)) return // no change
    demean = (demean0!=0)
    clear1()
}

transmorphic mm_qr::collin(| real scalar collin0)
{
    if (args()==0) return(collin)
    if (collin==(collin0!=0)) return // no change
    collin = (collin0!=0)
    clear2()
}

transmorphic mm_qr::p(| real scalar p0)
{
    if (args()==0) return(p)
    if (p0<=0 | p0>=1) _error(3300)
    if (p==p0) return // no change
    p = p0
    clear2b()
}

transmorphic mm_qr::b_init(| real colvector b_init0)
{
    if (args()==0) {
        if (b_init==.z) {
            err_nodata()
            _b_init()
        }
        return(b_init)
    }
    if (b_init0==.z) {  // use b_init(.z) to clear starting values
        b_init_user = 0
    }
    else {
        err_nodata()
        if (rows(b_init0)!=K) _error(3200) // wrong number of parameters
        b_init = b_init0
        b_init_user = 1
    }
    clear3()
}

transmorphic mm_qr::tol(| real scalar tol0)
{
    if (args()==0) return(tol)
    if (tol0<=0) _error(3300)
    if (tol==tol0) return // no change
    tol = tol0
    clear3()
}

transmorphic mm_qr::maxiter(| real scalar maxiter0)
{
    if (args()==0) return(maxiter)
    if (maxiter0<0) _error(3300)
    if (maxiter==maxiter0) return // no change
    maxiter = maxiter0
    clear3()
}

transmorphic mm_qr::beta(| real scalar beta0)
{
    if (args()==0) return(beta)
    if (beta0<0 | beta0>1) _error(3300)
    if (beta==beta0) return // no change
    beta = beta0
    clear3()
}

transmorphic mm_qr::method(| string scalar method0)
{
    if (args()==0) return(method)
    if (method==method0) return // no change
    if (!anyof(("fnb"), method0)) _error(3300)
    method = method0
    clear3()
}

transmorphic mm_qr::log(| real scalar log0)
{
    if (args()==0) return(log)
    if (!anyof((0,1,2,3), log0)) _error(3300)
    log = log0
}

real scalar mm_qr::n() return(n)

real scalar mm_qr::N() return(N)

real scalar mm_qr::cons() return(cons)

real scalar mm_qr::k() return(k)

real scalar mm_qr::K() return(K)

real colvector mm_qr::omit()
{
    if (omit==.z) {
        if (y!=NULL) set_collin()
    }
    return(omit)
}

real scalar mm_qr::k_omit()
{
    if (kadj==.z) {
        if (y!=NULL) set_collin()
    }
    return(K-Kadj)
}

real scalar mm_qr::ymean()
{
    if (ymean==.z) {
        if (y!=NULL) ymean = quadcross(*w, *y) / N
    }
    return(ymean)
}

real rowvector mm_qr::means()
{
    if (means==.z) {
        if (y!=NULL) means = quadcross(*w, *X) / N
    }
    return(means)
}

real colvector mm_qr::b()
{
    if (b==.z) {
        err_nodata()
        fit()
    }
    return(b)
}

real colvector mm_qr::xb(| real matrix X0)
{
    if (b==.z) (void) b()
    if (args()==1) {
        if (cols(X0)!=k) _error(3200)
        return(_xb(X0, b))
    }
    return(_xb(*X, b))
}

real colvector mm_qr::_xb(real matrix X, real colvector b)
{
    real scalar k
    
    k = rows(b) - cons
    if (k==0) {
        if (cons) return(J(rows(*y), 1, b))
        return(J(rows(*y), 1, .)) // model has no parameters
    }
    if (cons) return(X * b[|1\k|] :+ b[k+1])
    return(X * b)
}

real scalar mm_qr::gap()
{
    if (b==.z) (void) b()
    return(gap)
}

real scalar mm_qr::sdev()
{
    if (b==.z) (void) b()
    if (sdev==.z) sdev = _sdev(b)
    return(sdev)
}

real scalar mm_qr::_sdev(real colvector b)
{
    real colvector e
    
    e = *y - _xb(*X, b)
    return(cross(*w, (p :- (e:<0)) :* e))
}

real scalar mm_qr::iter()
{
    if (b==.z) (void) b()
    return(iter)
}

real scalar mm_qr::converged()
{
    if (b==.z) (void) b()
    return(conv)
}

void mm_qr::fit()
{
    real colvector b
    
    // identify collinear variables and exit if nothing to do
    set_collin()
    iter = conv = 0
    if (Kadj==0) { // model has no (non-omitted) parameters
        this.b = J(K,1,.)
        return
    }
    if (n==0) { // no observations
        this.b = J(K,1,.)
        return
    }
    // starting values
    if (b_init==.z) _b_init()
    b = -rmomit(meanadj(b_init, -1))
    
    // interior point algorithm (b will be replaced by solution)
    if      (method=="fnb")  fnb(b)
    else                     _error(3300)
    
    // rescale coefficients
    this.b = meanadj(addomit(-b), 1)
}

void mm_qr::set_collin()
{
    real scalar  nomit
    transmorphic t
    
    if (kadj!=.z) return    // already applied
    if (collin==0 | k==0) {
        kadj  = k
        Kadj  = K
        omit  = J(K, 1, 0)
        xindx = J(1, 0, .)
        return
    }
    t = lsfit()
    omit = mm_ls_omit(t)
    if (cons) omit[K] = 0 // not necessary, I believe
    nomit = sum(omit)
    kadj = k - nomit
    Kadj = K - nomit
    if (!nomit) xindx = J(1, 0, .)
    else        xindx = select(1..k, !omit[|1\k|]')
}

transmorphic mm_qr::lsfit()
{
    transmorphic t
    
    t = mm_ls(*y, *X, *w, cons, qd, demean)
    // hold on to coefficients for b_init
    b_ls = mm_ls_b(t)
    // hold on to means for demeaning
    if (cons & demean) { 
        ymean = mm_ls_ymean(t)
        means = mm_ls_means(t)
    }
    return(t)
}

void mm_qr::_b_init()
{
    if (K==0) { // model has no parameters
        b_init = J(K,1,.)
        return
    }
    if (n==0) { // no observations
        b_init = J(K,1,.)
        return
    }
    if (b_ls==.z) (void) lsfit()
    b_init = b_ls
}

real colvector mm_qr::rmomit(real colvector b)
{
    if (length(xindx)==0) return(b)
    return(b[xindx] \ (cons ? b[K] : J(0,1,.)))
}

real colvector mm_qr::addomit(real colvector b0)
{
    real colvector b
    
    if (length(xindx)==0) return(b0)
    b = J(K,1,0)
    b[xindx] = b0[|1\kadj|]
    if (cons) b[K] = b0[Kadj]
    return(b)
}

real colvector mm_qr::meanadj(real colvector b0, real scalar sign)
{
    real colvector b
    
    if (!demean) return(b0)
    if (!cons)   return(b0)
    b = b0
    b[K] = b[K] + ymean() * sign
    if (k) b[K] = b[K] - (means()*b[|1\k|]) * sign
    return(b)
}

// translation of rqfnb.f from "quantreg" package version 5.85 for R
void mm_qr::fnb(real colvector y)
{
    real scalar    fp, fd, mu, g
    real colvector b, dy, y0
    real colvector W, x, s, z, w, dx, ds, dz, dw, dxdz, dsdw, rhs, d, r
    real matrix    ada
    pointer scalar a
    
    // data
    s = -(*this.y)
    if (cons & demean)    s = s :+ ymean()
    a = get_X()
    if (rows(*this.w)!=1) W = (*this.w) * (n/N)
    else                  W = 1
    
    // algorithm
    gen_z_w(z=., w=., s - _xb(*a, y), tol)
    x = J(n, 1, 1-p)
    s = J(n, 1, p)
    if (maxiter<=0) {
        // make sure gap is filled in even if maxiter=0
        gap = cross(z, W, x) + cross(w, W, s)
    }
    b = cross(*a,cons, W, x,0)
    while (iter < maxiter) {
        iter = iter + 1
        d    = 1 :/ (z :/ x + w :/ s)
        ada  = cross(*a,cons, d:*W, *a,cons)
        r    = z - w
        rhs  = b - cross(*a,cons, W, x :- d:*r,0)
        dy   = cholsolve(ada, rhs)
        if (hasmissing(dy)) dy = invsym(ada) * rhs // singularity encountered
        dx   = d :* (_xb(*a, dy) - r)
        ds   = -dx
        dz   = -z :* (1 :+ dx :/ x)
        dw   = -w :* (1 :+ ds :/ s)
        fp   = minselect(x, dx, s, ds)
        fd   = minselect(w, dw, z, dz)
        if (min((fp, fd)) < 1) {
            mu   = cross(z, W, x) + cross(w, W, s)
            g    = cross(z + fd * dz, W, x + fp * dx) + 
                   cross(w + fd * dw, W, s + fp * ds)
            mu   = mu * (g / mu)^3 / (2 * n)
            dxdz = dx :* dz
            dsdw = ds :* dw
            r    = mu:/s - mu:/x + dxdz:/x - dsdw:/s
            rhs  = rhs + cross(*a,cons, d:*W, r,0)
            dy = cholsolve(ada, rhs)
            if (hasmissing(dy)) dy = invsym(ada) * rhs
            dx   = d :* (_xb(*a, dy) - z + w - r) 
            ds   = -dx
            dz   = (mu :- z :* dx :- dxdz) :/ x - z
            dw   = (mu :- w :* ds :- dsdw) :/ s - w
            fp   = minselect(x, dx, s, ds)
            fd   = minselect(w, dw, z, dz)
        }
        x    = x + fp * dx
        s    = s + fp * ds
        w    = w + fd * dw
        z    = z + fd * dz
        y0   = y
        y    = y + fd * dy
        gap = cross(z, W, x) + cross(w, W, s)
        if (log) {
            if (log>=2) printflush(".")
            else        printlog(mreldif(y, y0), iter)
        }
        if (gap<tol) {
            conv = 1
            break
        }
    }
    if (log==2 & iter) printflush("\n")
}

pointer scalar mm_qr::get_X()
{
    real matrix x
    
    if (!(cons & demean) & !length(xindx)) return(X)
    x = *X
    if (cons & demean) x = x :- means()
    if (length(xindx)) x = x[,xindx]
    return(&x)
}

void mm_qr::gen_z_w(real colvector z, real colvector w, real colvector r,
    real scalar eps)
{
    real colvector o
    
    z = r :* (r :> 0)
    w = z - r
    o = eps :* (abs(r):<eps)
    z = z + o
    w = w + o
}

real matrix mm_qr::cross(real matrix a, real matrix b, | real matrix c,
    real matrix d, real matrix e)
{
    if (args()==2) {
        if (qd) return(quadcross(a, b))
        return(::cross(a, b))
    }
    if (args()==3) {
        if (qd) return(quadcross(a, b, c))
        return(::cross(a, b, c))
    }
    if (args()==4) {
        if (qd) return(quadcross(a, b, c, d))
        return(::cross(a, b, c, d))
    }
    if (qd) return(quadcross(a, b, c, d, e))
    return(::cross(a, b, c, d, e))
}

real scalar mm_qr::minselect(real colvector x, real colvector dx, 
    real colvector s, real colvector ds)
{
    real colvector min, p
    
    min = J(3,1,.)
    p = select(1::n, dx:<0)
    if (length(p)) min[1] = min(-x[p] :/ dx[p]) * beta
    p = select(1::n, ds:<0)
    if (length(p)) min[2] = min(-s[p] :/ ds[p]) * beta
    min[3] = 1
    return(min(min))
}

// - wrapper for quick QR fit
real colvector mm_qrfit(
    real colvector      y,
    | real matrix       X, 
      real colvector    w, 
      real scalar       p,
      real scalar       cons,
      real colvector    b_init,
      real scalar       relax)
{
    class mm_qr scalar S
    
    S.data(y, X, w, cons)
    if (p!=.) S.p(p)
    if (args()>=6 & b_init!=.) S.b_init(b_init)
    if (args()<7) relax = 0
    if (!S.converged() & !relax) _error(3360)
    return(S.b())
}

end


*! {smcl}
*! {marker mm_mloc}{bf:mm_mloc.mata}{asis}
*! version 1.0.3  12mar2021  Ben Jann
version 11
mata:

// M-estimate of location
struct _mm_mls_struct {
    real scalar eff   // not used by mscale
    real scalar bp    // only used by mscale
    real scalar delta // only used by mscale
    real scalar k
    real scalar b
    real scalar b0
    real scalar s     // not used by mscale
    real scalar l     // only used by mscale
    real scalar d
    real scalar tol
    real scalar iter
    real scalar maxiter
    real scalar trace
    real scalar conv
}

struct _mm_mls_struct scalar mm_mloc(real colvector x,
  | real colvector w, real scalar eff, string scalar obj,
    real scalar b, real scalar s, 
    real scalar trace, real scalar tol, real scalar iter)
{
    real scalar b0, i
    real colvector z
    struct _mm_mls_struct scalar S
    pointer(real scalar function) scalar f
    
    if (args()<2) w = 1
    if (hasmissing(w)) _error(3351)
    if (hasmissing(x)) _error(3351)
    S.eff     = (eff<.   ? eff   : 95)
    S.trace   = (trace<. ? trace : 0)
    S.tol     = (tol<.   ? tol   : 1e-10)
    S.maxiter = (iter<.  ? iter  : st_numscalar("c(maxiter)"))
    S.iter    = 0
    S.conv    = 0
    S.b       = .
    if (obj=="" | obj=="huber") {
        f = &mm_huber_w()
        S.k = mm_huber_k(S.eff)
    }
    else if (obj=="biweight") {
        f = &mm_biweight_w()
        S.k = mm_biweight_k(S.eff)
    }
    else {
        printf("{err}'%s' not supported; invalid objective function\n", obj)
        _error(3498)
    }
    if (rows(x)==0) return(S) // nothing to do (b = ., conv = 0)
    S.b = S.b0 = (b<. ? b : mm_median(x, w))
    S.s  = (s<. ? s : mm_median(abs(x :- (b<. ? mm_median(x, w) : S.b)), w)
                      / invnormal(0.75))
    if (S.s<=0) return(S) // scale is 0 (b = starting value, conv = 0)
    // iterate
    z = (x :- S.b) / S.s
    for (i=1; i<=S.maxiter; i++) {
        b0 = S.b
        z = (*f)(z, S.k) :* w
        S.b = mean(x, z)
        S.d = abs(S.b - b0) / S.s
        if (S.trace) printf("{txt}{lalign 16:Iteration %g:}" +
            "absolute difference = %9.0g\n", i, S.d)
        if (S.d <= S.tol) break
        z = (x :- S.b) / S.s
    }
    S.iter = i
    S.conv = (S.iter<=S.maxiter)
    return(S)
}

real scalar    mm_mloc_b(struct _mm_mls_struct scalar S) return(S.b)
real scalar   mm_mloc_b0(struct _mm_mls_struct scalar S) return(S.b0)
real scalar    mm_mloc_s(struct _mm_mls_struct scalar S) return(S.s)
real scalar mm_mloc_conv(struct _mm_mls_struct scalar S) return(S.conv)
real scalar    mm_mloc_d(struct _mm_mls_struct scalar S) return(S.d)
real scalar mm_mloc_iter(struct _mm_mls_struct scalar S) return(S.iter)
real scalar    mm_mloc_k(struct _mm_mls_struct scalar S) return(S.k)
real scalar  mm_mloc_eff(struct _mm_mls_struct scalar S) return(S.eff)

// M-estimate of scale
struct _mm_mls_struct scalar mm_mscale(real colvector x,
  | real colvector w, real scalar bp,
    real scalar b, real scalar l, 
    real scalar trace, real scalar tol, real scalar iter)
{
    real scalar b0, i
    real colvector z
    struct _mm_mls_struct scalar S
    
    if (args()<2) w = 1
    if (hasmissing(w)) _error(3351)
    if (hasmissing(x)) _error(3351)
    S.bp      = (bp<.    ? bp    : 50)
    S.trace   = (trace<. ? trace : 0)
    S.tol     = (tol<.   ? tol   : 1e-10)
    S.maxiter = (iter<.  ? iter  : st_numscalar("c(maxiter)"))
    S.iter    = 0
    S.conv    = 0
    S.b       = .
    S.k       = mm_biweight_k_bp(S.bp)
    if (rows(x)==0) return(S) // nothing to do (b = ., conv = 0)
    S.l = (l<. ? l : mm_median(x, w))
    S.b = S.b0 = (b<. ? b : 
        mm_median(abs(x :- (l<. ? mm_median(x, w) : S.l)), w) / invnormal(0.75))
    if (S.b<=0) return(S) // scale is 0 (b = starting value, conv = 0)
    // iterate
    S.delta = S.bp/100 * S.k^2/6
    z = (x :- S.l)
    for (i=1; i<=S.maxiter; i++) {
        b0 = S.b
        S.b = sqrt(mean(mm_biweight_rho(z/b0, S.k), w) / S.delta) * b0
        S.d = abs(S.b/b0 - 1)
        if (S.trace) printf("{txt}{lalign 16:Iteration %g:}" +
            "absolute relative difference = %9.0g\n", i, S.d)
        if (S.d <= S.tol) break
    }
    S.iter = i
    S.conv = (S.iter<=S.maxiter)
    return(S)
}

real scalar     mm_mscale_b(struct _mm_mls_struct scalar S) return(S.b)
real scalar    mm_mscale_b0(struct _mm_mls_struct scalar S) return(S.b0)
real scalar     mm_mscale_l(struct _mm_mls_struct scalar S) return(S.l)
real scalar  mm_mscale_conv(struct _mm_mls_struct scalar S) return(S.conv)
real scalar     mm_mscale_d(struct _mm_mls_struct scalar S) return(S.d)
real scalar  mm_mscale_iter(struct _mm_mls_struct scalar S) return(S.iter)
real scalar     mm_mscale_k(struct _mm_mls_struct scalar S) return(S.k)
real scalar    mm_mscale_bp(struct _mm_mls_struct scalar S) return(S.bp)
real scalar mm_mscale_delta(struct _mm_mls_struct scalar S) return(S.delta)

// Huber tuning constant for given efficiency
real scalar mm_huber_k(real scalar eff)
{
    if (eff==95) return(1.34499751)
    if (eff==90) return( .98180232)
    if (eff==85) return( .73173882)
    if (eff==80) return( .52942958)
    if (eff<63.7)  _error(3498, "efficiency may not be smaller than 63.7")
    if (eff>99.9) _error(3498, "efficiency may not be larger than 99.9")
    return(round(mm_finvert(eff/100, &mm_huber_eff(), 0.001, 3), 1e-8))
}

// Huber efficiency for given tuning constant
real scalar mm_huber_eff(real scalar k)
{
    if (k<=0) return(2/pi())
    if (k<1e-4) {
        // use linear interpolation at bottom
        return(max((2/pi(), mm_huber_eff(1e-4) - (1e-4 - k)/1e-4 *
            (mm_huber_eff(2e-4)-mm_huber_eff(1e-4)))))
    }
    return((normal(k)-normal(-k))^2 / 
        (2 * (k^2 * (1 - normal(k)) + normal(k) - 0.5 - k * normalden(k))))
}

// Huber objective functions and weights
real colvector mm_huber_rho(real colvector x, real scalar k)
{
    real colvector y, d
    
    y = abs(x)
    d = y:<=k
    return((0.5 * x:^2):*d :+ (k*y :- 0.5*k^2):*(1 :- d))
}

real colvector mm_huber_psi(real colvector x, real scalar k)
{
    real colvector d
    
    d = abs(x):<=k
    return(x:*d :+ (sign(x)*k):*(1:-d))
}

real colvector mm_huber_phi(real colvector x, real scalar k)
{
    return(abs(x):<=k)
}

real colvector mm_huber_w(real colvector x, real scalar k)
{
    real colvector y
    
    y = abs(x)
    return(editmissing(k :/ y, 0):^(y:>k))
}

// biweight tuning constant for given efficiency
real scalar mm_biweight_k(real scalar eff0)
{
    real scalar k, eff
    
    if (eff0==95) return(4.6850649)
    if (eff0==90) return(3.8826616)
    if (eff0==85) return(3.4436898)
    if (eff0==80) return(3.1369087)
    if (eff0<0.1)  _error(3498, "efficiency may not be smaller than 0.1")
    if (eff0>99.9) _error(3498, "efficiency may not be larger than 99.9")
    eff = eff0/100
    k = 0.8376 + 1.499*eff + 0.7509*sin(1.301*eff^6) + 
        0.04945*eff/sin(3.136*eff) + 0.9212*eff/cos(1.301*eff^6)
    return(round(mm_finvert(eff, &mm_biweight_eff(), k/5, k*1.1), 1e-7))
}

// biweight efficiency for given tuning constant
real scalar mm_biweight_eff(real scalar k)
{   // using Simpson's rule integration
    real scalar    l, u, n, d, phi, psi2
    real colvector x, w
    
    if (k<=0)  return(0)
    if (k>100) return(1)
    l = 0; u = k; n = 1000; d = (u-l)/n 
    x = rangen(l, u, n+1)
    w = 1 \ colshape((J(n/2,1,4), J(n/2,1,2)), 1)
    w[n+1] = 1
    phi = 2 * (d / 3) * quadcolsum(_mm_biweight_eff_phi(x, k):*w)
    psi2 = 2 * (d / 3) * quadcolsum(_mm_biweight_eff_psi2(x, k):*w)
    return(phi^2 / psi2)
}

real matrix _mm_biweight_eff_phi(real matrix x, real scalar k)
{
    real matrix x2

    x2 = (x / k):^2
    return(normalden(x) :* ((1 :- x2) :* (1 :- 5*x2)) :* (x2:<=1))
}

real matrix _mm_biweight_eff_psi2(real matrix x, real scalar k)
{
    real matrix x2

    x2 = (x / k):^2
    return(normalden(x) :* ((x :* (1 :- x2):^2) :* (x2:<=1)):^2)
}

// biweight tuning constant for given breakdown point
real scalar mm_biweight_k_bp(real scalar bp)
{
    if (bp==50) return(1.547645)
    if (bp<1)  _error(3498, "bp may not be smaller than 1")
    if (bp>50) _error(3498, "bp may not be larger than 50")
    return(round(mm_finvert(bp/100, &mm_biweight_bp(), 1.5, 18), 1e-7))
}

// biweight breakdown point for given tuning constant
real scalar mm_biweight_bp(real scalar k)
{   // using Simpson's rule integration
    real scalar    l, u, n, d
    real colvector x, w

    if (k<=0) return(1)
    l = 0; u = k; n = 1000; d = (u-l)/n
    x = rangen(l, u, n+1)
    w = 1 \ colshape((J(n/2,1,4), J(n/2,1,2)), 1)
    w[n+1] = 1
    return(2 * (normal(-k) +  d/3 *
        quadcolsum(normalden(x):*mm_biweight_rho(x, k):*w) / (k^2/6)))
}

// biweight objective functions and weights
real colvector mm_biweight_rho(real colvector x, real scalar k)
{
    real colvector x2

    x2 = (x / k):^2
    return(k^2/6 * (1 :- (1 :- x2):^3):^(x2:<=1))
}

real colvector mm_biweight_psi(real colvector x, real scalar k)
{
    real colvector x2

    x2 = (x / k):^2
    return((x :* (1 :- x2):^2) :* (x2:<=1))
}

real colvector mm_biweight_phi(real colvector x, real scalar k)
{
    real colvector x2

    x2 = (x / k):^2
    return(((1 :- x2) :* (1 :- 5*x2)) :* (x2:<=1))
}
real colvector mm_biweight_w(real colvector x, real scalar k)
{
    real colvector x2

    x2 = (x / k):^2
    return(((1 :- x2):^2) :* (x2:<=1))
}

end


*! {smcl}
*! {marker mm_hl}{bf:mm_hl.mata}{asis}
*! version 1.0.0  21oct2020  Ben Jann
version 11
mata:

// robust statistics based on pairwise comparisons (code adapted from
// robstat.ado, version 1.0.3)
// - HL: location
// - Qn: scale
// - mc (medcouple): skewness

// HL estimator

real scalar mm_hl(real colvector X, | real colvector w, real scalar fw, 
    real scalar naive)
{
    real colvector p
    
    if (args()<2) w = 1
    if (args()<3) fw = 0
    if (args()<4) naive = 0
    if (hasmissing(w)) _error(3351)
    if (any(w:<0)) _error(3498, "negative weights not allowed")
    if (rows(X)==0) return(.)
    if (hasmissing(X)) _error(3351)
    if (naive) {
        if (rows(w)==1) return(_mm_hl_naive(X))
        if (fw) return(_mm_hl_naive_fw(X, w))
        return(_mm_hl_naive_w(X, w))
    }
    if (mm_issorted(X)) {
        if (rows(w)==1) return(_mm_hl(X))
        if (fw) return(_mm_hl_fw(X, w))
        return(_mm_hl_w(X, w))
    }
    if (rows(w)==1) return(_mm_hl(sort(X,1)))
    p = order((X,w), (1,2))
    if (fw) return(_mm_hl_fw(X[p], w[p]))
    return(_mm_hl_w(X[p], w[p]))
}

real scalar _mm_hl(real colvector x) // no weights
{
    // the trick of this algorithm is to consider only elements that are on the 
    // right of the main diagonal
    real scalar     i, j, m, k, n, nl, nr, nL, nR, trial
    real colvector  xx, /*ww,*/ l, r, L, R
    
    n = rows(x)                    // dimension of search matrix
    if (n==1) return(x)            // returning observed value if n=1
    xx      = /*ww =*/ J(n, 1, .)  // temp vector for matrix elements
    l = L   = (1::n):+1            // indices of left boundary (old and new)
    r = R   = J(n, 1, n)           // indices of right boundary (old and new)
    nl = nl = n + comb(n, 2)       // number of cells below left boundary
    nr = nR = n * n                // number of cells within right boundary
    k       = nl + comb(n, 2)/2    // target quantile
    while ((nr-nl)>n) {
        // get trial value
        m = 0
        for (i=1; i<n; i++) { // last row cannot contain candidates
            if (l[i]<=r[i]) {
                // high median within row
                xx[++m] = __mm_hl_el(x, i, l[i]+trunc((r[i]-l[i]+1)/2))
                /*m++
                ww[m] = r[i] - l[i] + 1
                xx[m] = __mm_hl_el(x, i, l[i]+trunc(ww[m]/2))*/
            }
        }
        trial = _mm_hl_qhi(xx[|1 \ m|], .5)
        /*trial = _mm_hl_qhi_w(xx[|1 \ m|], ww[|1 \ m|], .5)*/
        /*the unweighted quantile is faster; results are the same*/
        // move right border
        j = n-1
        for (i=(n-1); i>=1; i--) {
            if (i==j) {
                if (__mm_hl_el(x, i, j)>=trial) {
                    R[i] = j
                    j = i-1
                    continue
                }
            }
            if (j<n) {
                while (__mm_hl_el(x, i, j+1)<trial) {
                    j++
                    if (j==n) break
                }
            }
            R[i] = j
        }
        nR = sum(R)
        if (nR>k) {
            swap(r, R)
            nr = nR
            continue
        }
        // move left border
        j = n + 1
        for (i=1; i<=n; i++) {
            if (j>(i+1)) {
                while (__mm_hl_el(x, i, j-1)>trial) {
                    j--
                    if (j==(i+1)) break
                }
            }
            if (j<(i+1)) j = (i+1)
            L[i] = j
        }
        nL = sum(L) - n
        if (nL<k) {
            swap(l, L)
            nl = nL
            continue
        }
        // trial = low quantile = high quantile
        if (ceil(k)!=k | (nR<k & nL>k)) return(trial)
        // trial = low quantile
        if (nL==k) {
            m = 0
            for (i=1; i<=n; i++) {
                if (L[i]>n) continue
                xx[++m] = __mm_hl_el(x, i, L[i])
            }
            return((trial+min(xx[|1 \ m|]))/2)
        }
        // trial = high quantile
        m = 0
        for (i=1; i<=n; i++) {
            if (R[i]<=i) continue
            xx[++m] = __mm_hl_el(x, i, R[i])
        }
        return((trial+max(xx[|1 \ m|]))/2)
    }
    // get target value from remaining candidates
    m = 0
    for (i=1; i<n; i++) { // last row cannot contain candidates
        if (l[i]<=r[i]) {
            for (j=l[i]; j<=r[i]; j++) {
                m++
                xx[m] = __mm_hl_el(x, i, j)
            }
        }
    }
    return(_mm_hl_q(xx[|1 \ m|], k, nl))
}

real scalar __mm_hl_el(real colvector y, real scalar i, real scalar j)
{
    return((y[i] + y[j])/2)
}

real scalar _mm_hl_w(real colvector x, real colvector w)
{
    real scalar     i, j, m, k, n, nl, nr, trial, Wl, WR, WL, W0, W1
    real colvector  xx, ww, l, r, L, R, ccw
    
    n       = rows(x)              // dimension of search matrix
    if (n==1) return(x)            // returning observed value if n=1
    xx = ww = J(n, 1, .)           // temp vector for matrix elements
    l = L   = (1::n):+1            // indices of left boundary (old and new)
    r = R   = J(n, 1, n)           // indices of right boundary (old and new)
    nl      = comb(n, 2) + n       // number of cells below left boundary
    nr      = n * n                // number of cells within right boundary
    ccw     = quadrunningsum(w)    // cumulative column weights
    W1      = quadsum(w[|2 \ .|] :* ccw[|1 \ n-1|]) // sum of weights in target triangle
    W0      = W1 + quadsum(w:*w)   // sum of weights in rest of search matrix
    Wl = WL = W0                   // sum of weights below left boundary
    WR      = W0 + W1              // sum of weights within right boundary
    k       = W0 + W1/2            // target quantile
    while ((nr-nl)>n) {
        // get trial value
        m = 0
        for (i=1; i<n; i++) { // last row cannot contain candidates
            if (l[i]<=r[i]) {
                // high median within row
                xx[++m] = __mm_hl_el(x, i, l[i]+trunc((r[i]-l[i]+1)/2))
            }
        }
        trial = _mm_hl_qhi(xx[|1 \ m|], .5)
        // move right border
        j = n-1
        for (i=(n-1); i>=1; i--) {
            if (i==j) {
                if (__mm_hl_el(x, i, j)>=trial) {
                    R[i] = j
                    j = i-1
                    continue
                }
            }
            if (j<n) {
                while (__mm_hl_el(x, i, j+1)<trial) {
                    j++
                    if (j==n) break
                }
            }
            R[i] = j
        }
        WR = quadsum(w:*ccw[R])
        if (WR>k) {
            swap(r, R)
            nr = sum(R)
            continue
        }
        // move left border
        j = n + 1
        for (i=1; i<=n; i++) {
            if (j>(i+1)) {
                while (__mm_hl_el(x, i, j-1)>trial) {
                    j--
                    if (j==(i+1)) break
                }
            }
            if (j<(i+1)) j = (i+1)
            L[i] = j
        }
        WL = quadsum(w:*ccw[L:-1])
        if (WL<k) {
            swap(l, L)
            Wl = WL
            nl = sum(L) - n
            continue
        }
        // trial = low quantile = high quantile
        if (WR==WL | (WR<k & WL>k)) return(trial)
        // trial = low quantile
        if (WL==k) {
            m = 0
            for (i=1; i<=n; i++) {
                if (L[i]>n) continue
                xx[++m] = __mm_hl_el(x, i, L[i])
            }
            return((trial+min(xx[|1 \ m|]))/2)
        }
        // trial = high quantile
        m = 0
        for (i=1; i<=n; i++) {
            if (R[i]<=i) continue
            xx[++m] = __mm_hl_el(x, i, R[i])
        }
        return((trial+max(xx[|1 \ m|]))/2)
    }
    // get target value from remaining candidates
    m = 0
    for (i=1; i<n; i++) { // last row cannot contain candidates
        if (l[i]<=r[i]) {
            for (j=l[i]; j<=r[i]; j++) {
                m++
                xx[m] = __mm_hl_el(x, i, j)
                ww[m] = w[i] * w[j]
            }
        }
    }
    return(_mm_hl_q_w(xx[|1 \ m|], ww[|1 \ m|], k, Wl))
}

real scalar _mm_hl_fw(real colvector x, real colvector w)
{   // the algorithm "duplicates" the diagonal so that the relevant pairs in 
    // case of w>1 can be taken into account
    real scalar     i, j, m, k, n, nl, nr, trial, Wl, WR, WL, W0, W1
    real colvector  xx, ww, l, r, L, R, ccw, wcorr, idx
    
    if (any(trunc(w):!=w)) _error(3498, "non-integer frequency not allowed")
    n       = rows(x)          // dimension of search matrix
    if (n==1) return(x)        // returning observed value if n=1
    xx = ww = J(n, 1, .)       // temp vector for matrix elements
    ccw     = runningsum(w)    // cumulative column weights
    wcorr   = mm_cond(w:<=1, 0, comb(w, 2)) // correction of weights
    idx     = 1::n             // diagonal indices
    l = L   = idx :+ 1 :+ (wcorr:==0) // indices of left boundary (old and new)
    r = R   = J(n, 1, n+1)     // indices of right boundary (old and new)
    nl      = comb(n, 2) + n + sum(wcorr:==0) // n. of cells below left boundary
    nr      = n * (n+1)        // number of cells within right boundary
    W0 = W1 = sum(w[|2 \ .|] :* ccw[|1 \ n-1|])
    W1      = W1 + sum(wcorr)  // sum of weights in target triangle
    W0      = W0 + sum(w) + sum(wcorr) // sum of weights in rest of search matrix
    Wl = WL = W0               // sum of weights below left boundary
    WR      = W0 + W1          // sum of weights within right boundary
    k       = W0 + W1/2        // target quantile
    while ((nr-nl)>n) {
        // get trial value
        m = 0
        for (i=1; i<=n; i++) {
            if (l[i]<=r[i]) {
                // high median within row
                xx[++m] = __mm_hl_el(x, i, 
                    (l[i]-1) + trunc(((r[i]-1)-(l[i]-1)+1)/2))
            }
        }
        trial = _mm_hl_qhi(xx[|1 \ m|], .5)
        // move right border
        j = n
        for (i=n; i>=1; i--) {
            if (j==i) {
                if (__mm_hl_el(x, i, j)>=trial) {
                    R[i] = j + (wcorr[i]==0)
                    j = i-1
                    continue
                }
            }
            if (j<=n) {
                while (__mm_hl_el(x, i, (j+1)-((j+1)>i))<trial) {
                    j++
                    if (j>n) break
                }
            }
            R[i] = j
        }
        WR = sum(w:*ccw[R:-(R:>idx)]) - sum(wcorr:*(R:==idx))
        if (WR>k) {
            swap(r, R)
            nr = sum(R)
            continue
        }
        // move left border
        j = n + 2
        for (i=1; i<=n; i++) {
            if (j>(i+1)) {
                while (__mm_hl_el(x, i, (j-1)-((j-1)>i))>trial) {
                    j--
                    if (j==(i+1)) break
                }
            }
            if (j<(i+1)) j = (i+1)
            if (j==(i+1)) {
                if (wcorr[i]==0) j++
            }
            L[i] = j
        }
        WL = sum(w:*ccw[L:-1:-((L:-1):>idx)]) - sum(wcorr:*((L:-1):==idx))
        if (WL<k) {
            swap(l, L)
            Wl = WL
            nl = sum(L) - n
            continue
        }
        // trial = low quantile = high quantile
        if (ceil(k)!=k | (WR<k & WL>k)) return(trial)
        // trial = low quantile
        if (WL==k) {
            m = 0
            for (i=1; i<=n; i++) {
                if ((L[i]-(L[i]>i))>n) continue
                xx[++m] = __mm_hl_el(x, i, L[i]-(L[i]>i))
            }
            return((trial+min(xx[|1 \ m|]))/2)
        }
        // trial = high quantile
        m = 0
        for (i=1; i<=n; i++) {
            if (R[i]<=i) continue
            if (R[i]==i+1) {
                if (wcorr[i]==0) continue
            }
            xx[++m] = __mm_hl_el(x, i, R[i]-1)
        }
        return((trial+max(xx[|1 \ m|]))/2)
    }
    // get target value from remaining candidates
    m = 0
    for (i=1; i<=n; i++) {
        if (l[i]<=r[i]) {
            for (j=l[i]; j<=r[i]; j++) {
                m++
                xx[m] = __mm_hl_el(x, i, j-1)
                ww[m] = (i==(j-1) ? comb(w[i], 2) : w[i]*w[j-1])
            }
        }
    }
    return(_mm_hl_q_w(xx[|1 \ m|], ww[|1 \ m|], k, Wl))
}

real scalar _mm_hl_naive(real colvector x) // no weights
{
    real scalar    i, j, m, n
    real colvector xx
    
    n = rows(x)
    if (n==1) return(x) // HL undefined if n=1; returning observed value
    m = 0
    xx = J(comb(n,2), 1, .)
    for (i=1; i<n; i++) {
        for (j=(i+1); j<=n; j++) {
            xx[++m] = __mm_hl_el(x, i, j)
        }
    }
    return(_mm_hl_q(xx, .5))
}

real scalar _mm_hl_naive_w(real colvector x, real colvector w)
{
    real scalar    i, j, m, n
    real colvector xx, ww

    n = rows(x)
    if (n==1) return(x) // HL undefined if n=1; returning observed value
    m = 0
    xx = J(comb(n,2), 1, .)
    ww = J(rows(xx), 1, .)
    for (i=1; i<n; i++) {
        for (j=(i+1); j<=n; j++) {
            m++
            xx[m] = __mm_hl_el(x, i, j)
            ww[m] = w[i]*w[j]
        }
    }
    return(_mm_hl_q_w(xx, ww, .5))
}

real scalar _mm_hl_naive_fw(real colvector x, real colvector w)
{
    real scalar    i, j, m, n
    real colvector xx, ww

    if (any(trunc(w):!=w)) _error(3498, "non-integer frequency not allowed")
    n = rows(x)
    if (n==1) return(x) // HL undefined if n=1; returning observed value
    m = 0
    xx = J(comb(n,2)+sum(w:>1), 1, .)
    ww = J(rows(xx), 1, .)
    for (i=1; i<=n; i++) {
        for (j=i; j<=n; j++) {
            if (i==j) {
                if (w[i]==1) continue
            }
            m++
            xx[m] = __mm_hl_el(x, i, j)
            ww[m] = (i==j ? comb(w[i], 2) : w[i]*w[j])
        }
    }
    return(_mm_hl_q_w(xx, ww, .5))
}

// Qn estimator

real scalar mm_qn(real colvector X, | real colvector w, real scalar fw, 
    real scalar naive)
{
    real colvector p
    real scalar    c
    
    if (args()<2) w = 1
    if (args()<3) fw = 0
    if (args()<4) naive = 0
    if (hasmissing(w)) _error(3351)
    if (any(w:<0)) _error(3498, "negative weights not allowed")
    if (rows(X)==0) return(.)
    if (hasmissing(X)) _error(3351)
    c = 1 / (sqrt(2) * invnormal(5/8))
    if (naive) {
        if (rows(w)==1) return(_mm_qn_naive(X) * c)
        if (fw) return(_mm_qn_naive_fw(X, w) * c)
        return(_mm_qn_naive_w(X, w) * c)
    }
    if (mm_issorted(X)) {
        if (rows(w)==1) return(_mm_qn(X) * c)
        if (fw) return(_mm_qn_fw(X, w) * c)
        return(_mm_qn_w(X, w) * c)
    }
    if (rows(w)==1) return(_mm_qn(sort(X,1)) * c)
    p = order((X,w), (1,2))
    if (fw) return(_mm_qn_fw(X[p], w[p]) * c)
    return(_mm_qn_w(X[p], w[p]) * c)
}

real scalar _mm_qn(real colvector x) // no weights
{
    real scalar     i, j, m, k, n, nl, nr, nL, nR, trial
    real colvector  xx, /*ww,*/ l, r, L, R

    n = rows(x)                    // dimension of search matrix
    if (n==1) return(0)            // returning zero if n=1
    xx      = /*ww =*/ J(n, 1, .)  // temp vector for matrix elements
    l = L   = (n::1):+1            // indices of left boundary (old and new)
    r = R   = J(n, 1, n)           // indices of right boundary (old and new)
    nl = nl = comb(n, 2) + n       // number of cells below left boundary
    nr = nR = n * n                // number of cells within right boundary
    k       = nl + comb(n, 2)/4    // target quantile
    /*k = nl + comb(trunc(n/2) + 1, 2)*/
    while ((nr-nl)>n) {
        // get trial value
        m = 0
        for (i=2; i<=n; i++) { // first row cannot contain candidates
            if (l[i]<=r[i]) {
                // high median within row
                xx[++m] = __mm_qn_el(x, i, l[i]+trunc((r[i]-l[i]+1)/2), n)
                /*m++
                ww[m] = r[i] - l[i] + 1
                xx[m] = __mm_qn_el(x, i, l[i]+trunc(www[m]/2), n)*/
            }
        }
        trial = _mm_hl_qhi(xx[|1 \ m|], .5)
        /*trial = _mm_hl_qhi_w(xx[|1 \ m|], ww[|1 \ m|], .5)*/
        /*the unweighted quantile is faster; results are the same*/
        //move right border
        j = 0
        for (i=n; i>=1; i--) {
            if (j<n) {
                while (__mm_qn_el(x, i, j+1, n)<trial) {
                    j++
                    if (j==n) break
                }
            }
            R[i] = j
        }
        nR = sum(R)
        if (nR>k) {
            swap(r, R)
            nr = nR
            continue
        }
        // move left border
        j = n + 1
        for (i=1; i<=n; i++) {
            while (__mm_qn_el(x, i, j-1, n)>trial) {
                j--
            }
            L[i] = j
        }
        nL = sum(L) - n
        if (nL<k) {
            swap(l, L)
            nl = nL
            continue
        }
        // trial = low quantile = high quantile
        if (ceil(k)!=k | (nR<k & nL>k)) return(trial)
        // trial = low quantile
        if (nL==k) {
            m = 0
            for (i=1; i<=n; i++) {
                if (L[i]>n) continue
                xx[++m] = __mm_qn_el(x, i, L[i], n)
            }
            return((trial+min(xx[|1 \ m|]))/2)
        }
        // trial = high quantile
        for (i=1; i<=n; i++) {
            xx[i] = __mm_qn_el(x, i, R[i], n)
        }
        return((trial+max(xx))/2)
    }
    // get target value from remaining candidates
    m = 0
    for (i=2; i<=n; i++) { // first row cannot contain candidates
        if (l[i]<=r[i]) {
            for (j=l[i]; j<=r[i]; j++) {
                m++
                xx[m] = __mm_qn_el(x, i, j, n)
            }
        }
    }
    return(_mm_hl_q(xx[|1 \ m|], k, nl))
}

real scalar __mm_qn_el(real colvector y, real scalar i, real scalar j, 
    real scalar n)
{
    return(y[i] - y[n-j+1])
}

real scalar _mm_qn_w(real colvector x, real colvector w)
{
    real scalar     i, j, m, k, n, nl, nr, trial, Wl, WR, WL, W0, W1
    real colvector  xx, ww, l, r, L, R, p, ccw

    n       = rows(x)              // dimension of search matrix
    if (n==1) return(0)            // returning zero if n=1
    xx = ww = J(n, 1, .)           // temp vector for matrix elements
    l = L   = (n::1):+1            // indices of left boundary (old and new)
    r = R   = J(n, 1, n)           // indices of right boundary (old and new)
    nl      = comb(n, 2) + n       // number of cells below left boundary
    nr      = n * n                // number of cells within right boundary
    ccw     = quadrunningsum(w[n::1]) // cumulative column weights
    W0      = quadsum(w:*ccw[n::1]) // sum weights in rest of search matrix
    W1      = quadsum(w:*ccw[rows(ccw)]) - W0 // sum of weights target triangle
    Wl = WL = W0                   // sum of weights below left boundary
    WR      = W0 + W1              // sum of weights within right boundary
    k       = W0 + W1/4            // target sum (high 25% quantile)
    while ((nr-nl)>n) {
        // get trial value
        m = 0
        for (i=2; i<=n; i++) { // first row cannot contain candidates
            if (l[i]<=r[i]) {
                // high median within row
                xx[++m] = __mm_qn_el(x, i, l[i]+trunc((r[i]-l[i]+1)/2), n)
            }
        }
        trial = _mm_hl_qhi(xx[|1 \ m|], .5)
        //move right border
        j = 0
        for (i=n; i>=1; i--) {
            if (j<n) {
                while (__mm_qn_el(x, i, j+1, n)<trial) {
                    j++
                    if (j==n) break
                }
            }
            R[i] = j
        }
        p = (R:>0)
        WR = quadsum(select(w, p) :* ccw[select(R, p)])
        if (WR>k) {
            swap(r, R)
            nr = sum(R)
            continue
        }
        // move left border
        j = n + 1
        for (i=1; i<=n; i++) {
            while (__mm_qn_el(x, i, j-1, n)>trial) {
                j--
            }
            L[i] = j
        }
        WL = quadsum(w :* ccw[L:-1])
        if (WL<k) {
            swap(l, L)
            Wl = WL
            nl = sum(L) - n
            continue
        }
        // trial = low quantile = high quantile
        if (WR==WL | (WR<k & WL>k)) return(trial)
        // trial = low quantile
        if (WL==k) {
            m = 0
            for (i=1; i<=n; i++) {
                if (L[i]>n) continue
                xx[++m] =  __mm_qn_el(x, i, L[i], n)
            }
            return((trial+min(xx[|1 \ m|]))/2)
        }
        // trial = high quantile
        for (i=1; i<=n; i++) {
            xx[i] = __mm_qn_el(x, i, R[i], n)
        }
        return((trial+max(xx))/2)
    }
    // get target value from remaining candidates
    m = 0
    for (i=2; i<=n; i++) { // first row cannot contain candidates
        if (l[i]<=r[i]) {
            for (j=l[i]; j<=r[i]; j++) {
                m++
                xx[m] = __mm_qn_el(x, i, j, n)
                ww[m] = w[i] * w[n-j+1]
            }
        }
    }
    return(_mm_hl_q_w(xx[|1 \ m|], ww[|1 \ m|], k, Wl))
}

real scalar _mm_qn_fw(real colvector x, real colvector w)
{
    real scalar     i, j, m, k, n, nl, nr, trial, Wl, WR, WL, W0, W1
    real colvector  xx, ww, l, r, L, R, p, ccw, wcorr, idx

    if (any(trunc(w):!=w)) _error(3498, "non-integer frequency not allowed")
    n       = rows(x)              // dimension of search matrix
    if (n==1) return(0)            // returning zero if n=1
    xx = ww = J(n, 1, .)           // temp vector for matrix elements
    idx     = n::1                 // (minor) diagonal indices
    l = L   = idx:+1               // indices of left boundary (old and new)
    r = R   = J(n, 1, n+1)         // indices of right boundary (old and new)
    nl      = comb(n, 2) + n       // number of cells below left boundary
    nr      = n * (n+1)            // number of cells within right boundary
    ccw     = runningsum(w[idx])   // cumulative column weights
    wcorr   = mm_cond(w:<=1, 0, comb(w, 2))[idx] // correction of weights
    W0      = sum(w:*ccw[idx]) - sum(wcorr) // sum of weights in rest of search matrix
    W1      = sum(w:*ccw[rows(ccw)]) - W0 // sum of weights in target triangle
    Wl = WL = W0                   // sum of weights below left boundary
    WR      = W0 + W1              // sum of weights within right boundary
    k       = W0 + W1/4            // target quantile
    while ((nr-nl)>n) {
        // get trial value
        m = 0
        for (i=1; i<=n; i++) {
            if (l[i]<=r[i]) {
                // high median within row
                xx[++m] = __mm_qn_el(x, i, 
                    (l[i]-1)+trunc(((r[i]-1)-(l[i]-1)+1)/2), n)
            }
        }
        trial = _mm_hl_qhi(xx[|1 \ m|], .5)
        //move right border
        j = 0
        for (i=n; i>=1; i--) {
            if (j<=n) {
                while (__mm_qn_el(x, i, (j+1)-((j+1)>(n-i+1)), n)<trial) {
                    j++
                    if (j>n) break
                }
            }
            R[i] = j
        }
        p = (R:>0)
        WR = sum(select(w, p) :* ccw[select(R:-(R:>idx), p)]) - sum(wcorr:*(R:==idx))
        if (WR>k) {
            swap(r, R)
            nr = sum(R)
            continue
        }
        // move left border
        j = n + 2
        for (i=1; i<=n; i++) {
            while (__mm_qn_el(x, i, (j-1)-((j-1)>(n-i+1)), n)>trial) {
                j--
            }
            L[i] = j
        }
        WL = sum(w:*ccw[L:-1:-((L:-1):>idx)]) - sum(wcorr:*((L:-1):==idx))
        if (WL<k) {
            swap(l, L)
            Wl = WL
            nl = sum(L) - n
            continue
        }
        // trial = low quantile = high quantile
        if (ceil(k)!=k | (WR<k & WL>k)) return(trial)
        // trial = low quantile
        if (WL==k) {
            m = 0
            for (i=1; i<=n; i++) {
                if ((L[i]-(L[i]>(n-i+1)))>n) continue
                xx[++m] = __mm_qn_el(x, i, L[i]-(L[i]>(n-i+1)), n)
            }
            return((trial+min(xx[|1 \ m|]))/2)
        }
        // trial = high quantile
        for (i=1; i<=n; i++) {
            xx[i] = __mm_qn_el(x, i, R[i]-(R[i]>(n-i+1)), n)
        }
        return((trial+max(xx))/2)
    }
    // get target value from remaining candidates
    m = 0
    for (i=1; i<=n; i++) {
        if (l[i]<=r[i]) {
            for (j=l[i]; j<=r[i]; j++) {
                m++
                xx[m] = __mm_qn_el(x, i, j-1, n)
                ww[m] = ((n-i+1)==(j-1) ? comb(w[i], 2) : w[i]*w[n-(j-1)+1])
            }
        }
    }
    return(_mm_hl_q_w(xx[|1 \ m|], ww[|1 \ m|], k, Wl))
}

real scalar _mm_qn_naive(real colvector x) // no weights
{
    real scalar    i, j, m, n
    real colvector xx
    
    n = rows(x)
    if (n==1) return(0) // returning zero if n=1
    m = 0
    xx = J(comb(n,2), 1, .)
    for (i=1; i<n; i++) {
        for (j=(i+1); j<=n; j++) {
            xx[++m] = abs(x[i] - x[j])
        }
    }
    return(_mm_hl_q(xx, 0.25))
}

real scalar _mm_qn_naive_w(real colvector x, real colvector w)
{
    real scalar    i, j, m, n
    real colvector xx, ww
    
    n = rows(x)
    if (n==1) return(0) // Qn undefined if n=1; returning zero
    m = 0
    xx = J(comb(n,2), 1, .)
    ww = J(rows(xx), 1, .)
    for (i=1; i<n; i++) {
        for (j=(i+1); j<=n; j++) {
            m++
            xx[m] = abs(x[i] - x[j])
            ww[m] = w[i]*w[j]
        }
    }
    return(_mm_hl_q_w(xx, ww, 0.25))
}

real scalar _mm_qn_naive_fw(real colvector x, real colvector w)
{
    real scalar    i, j, m, n
    real colvector xx, ww

    if (any(trunc(w):!=w)) _error(3498, "non-integer frequency not allowed")
    n = rows(x)
    if (n==1) return(0) // returning zero if n=1
    m = 0
    xx = J(comb(n,2)+sum(w:>1), 1, .)
    ww = J(rows(xx), 1, .)
    for (i=1; i<=n; i++) {
        for (j=i; j<=n; j++) {
            if (i==j) {
                if (w[i]==1) continue
            }
            m++
            xx[m] = abs(x[i] - x[j])
            ww[m] = (i==j ? comb(w[i], 2) : w[i]*w[j])
        }
    }
    return(_mm_hl_q_w(xx, ww, 0.25))
}

// medcouple

real scalar mm_mc(real colvector X, | real colvector w, real scalar fw, 
    real scalar naive)
{
    real colvector p, r
    
    if (args()<2) w = 1
    if (args()<3) fw = 0
    if (args()<4) naive = 0
    if (hasmissing(w)) _error(3351)
    if (any(w:<0)) _error(3498, "negative weights not allowed")
    if (rows(X)==0) return(.)
    if (hasmissing(X)) _error(3351)
    if (mm_issorted(X)) {
        if (naive) {
            if (rows(w)==1) return(_mm_mc_naive(X:-_mm_median(X)))
            if (fw) return(_mm_mc_naive_fw(X:-_mm_median(X,w), w))
            return(_mm_mc_naive_w(X:-_mm_median(X,w), w))
        }
        r = rows(X)::1
        if (rows(w)==1) return(_mm_mc(X[r]:-_mm_median(X)))
        if (fw) return(_mm_mc_fw(X[r]:-_mm_median(X,w), w[r]))
        return(_mm_mc_w(X[r]:-_mm_median(X,w), w[r]))
    }
    if (naive) {
        if (rows(w)==1) return(_mm_mc_naive(X:-mm_median(X)))
        if (fw) return(_mm_mc_naive_fw(X:-mm_median(X,w), w))
        return(_mm_mc_naive_w(X:-mm_median(X,w), w))
    }
    if (rows(w)==1) {
        p = order(X, 1)
        r = p[rows(X)::1]
        return(_mm_mc(X[r]:-_mm_median(X[p])))
    }
    p = order((X,w), (1,2))
    r = p[rows(X)::1]
    if (fw) return(_mm_mc_fw(X[r]:-_mm_median(X[p], w[p]), w[r]))
    return(_mm_mc_w(X[r]:-_mm_median(X[p], w[p]), w[r]))
}

real scalar _mm_mc(real colvector x) // no weights; assumes med(x)=0
{
    real scalar     i, j, m, k, n, q, nl, nr, nL, nR, trial, npos, nzero
    real colvector  xx, /*ww,*/ l, r, L, R, xpos, xneg

    if (rows(x)<=1) return(0)      // returning zero if n=1
    xpos    = select(x, x:>0)      // obervations > median
    npos    = rows(xpos)           // number of obs > median
    xneg    = select(x, x:<0)      // observations < median
    nzero   = sum(x:==0)           // number of obs = median
    n       = npos + nzero         // number of rows in search matrix
    q       = rows(xneg) + nzero   // number of columns in search matrix
    xx      = /*ww =*/ J(n, 1, .)  // temp vector for matrix elements
    l = L   = J(n, 1, 1)           // indices of left boundary (old and new)
    r = R   = J(n, 1, q)           // indices of right boundary (old and new)
    nl = nL = 0                    // number of cells below left boundary
    nr = nR = n * q                // number of cells within right boundary
    k       = n*q/2                // target quantile
    while ((nr-nl)>n) {
        // get trial value
        m = 0
        for (i=1; i<=n; i++) {
            if (l[i]<=r[i]) {
                // high median within row
                xx[++m] = -__mm_mc_el(xpos, xneg, npos, nzero, i,
                           l[i]+trunc((r[i]-l[i]+1)/2))
                /*m++
                ww[m] = r[i] - l[i] + 1
                xx[m] = -__mm_mc_el(xpos, xneg, npos, nzero, i,
                          l[i]+trunc((ww[m])/2))*/
            }
        }
        trial = _mm_hl_qhi(xx[|1 \ m|], .5)
        /*the unweighted quantile is faster; results are the same*/
        /*trial = _mm_hl_qhi_w(xx[|1 \ m|], ww[|1 \ m|], .5)*/
        // move right border
        j = 0
        for (i=n; i>=1; i--) {
            if (j<q) {
                while (-__mm_mc_el(xpos, xneg, npos, nzero, i, j+1)<trial) {
                    j++
                    if (j==q) break
                }
            }
            R[i] = j
        }
        nR = sum(R)
        if (nR>k) {
            swap(r, R)
            nr = nR
            continue
        }
        // move left border
        j = q + 1
        for (i=1; i<=n; i++) {
            while (-__mm_mc_el(xpos, xneg, npos, nzero, i, j-1)>trial) {
                j--
            }
            L[i] = j
        }
        nL = sum(L) - n
        if (nL<k) {
            swap(l, L)
            nl = nL
            continue
        }
        // trial = low quantile = high quantile
        if (ceil(k)!=k | (nR<k & nL>k)) return(-trial)
        // trial = low quantile
        if (nL==k) {
            m = 0
            for (i=1; i<=n; i++) {
                if (L[i]>q) continue
                xx[++m] = -__mm_mc_el(xpos, xneg, npos, nzero, i, L[i])
            }
            return(-(trial+min(xx[|1 \ m|]))/2)
        }
        // trial = high quantile
        for (i=1; i<=n; i++) {
            xx[i] = -__mm_mc_el(xpos, xneg, npos, nzero, i, R[i])
        }
        return(-(trial+max(xx))/2)
    }
    // get target value from remaining candidates
    m = 0
    for (i=1; i<=n; i++) {
        if (l[i]<=r[i]) {
            for (j=l[i]; j<=r[i]; j++) {
                m++
                xx[m] = -__mm_mc_el(xpos, xneg, npos, nzero, i, j)
            }
        }
    }
    return(-_mm_hl_q(xx[|1 \ m|], k, nl))
}

real scalar __mm_mc_el(real colvector xpos, real colvector xneg, 
    real scalar npos, real scalar nzero, real scalar i, real scalar j)
{
    if (i<=npos) {
        if (j<=nzero) return(1)
        // => (j>nzero)
        return((xpos[i] + xneg[j-nzero])/(xpos[i] - xneg[j-nzero]))
    }
    // => (i>npos)
    if (j>nzero)  return(-1)
    // => (j<=nzero)
    return(sign((npos+nzero-i+1)-j))
}

real scalar _mm_mc_w(real colvector x, real colvector w) // assumes med(x)=0
{
    real scalar     i, j, m, k, n, q, nl, nr, trial, npos, nzero, Wl, WR, WL, W
    real colvector  xx, ww, l, r, L, R, p, ccw, xpos, xneg, wpos, wneg, wzero

    if (rows(x)<=1) return(0)      // returning zero if n=1
    p = (x:>0)
    xpos    = select(x, p)         // obervations > median
    wpos    = select(w, p)         // weights of obervations >= median
    npos    = rows(xpos)           // number of obs > median
    p = (x:<0)
    xneg    = select(x, p)         // observations < median
    wneg    = select(w, p)         // weights of obervations <= median
    wzero   = select(w, x:==0)     // weights of obervations = median
    nzero   = rows(wzero)          // number of obs = median
    if (nzero>0) {
        wpos = wpos \ wzero
        wneg = wzero[nzero::1] \ wneg // need to use reverse ordered wzero
    }
    n       = npos + nzero         // number of rows in search matrix
    q       = rows(xneg) + nzero   // number of columns in search matrix
    xx = ww = J(n, 1, .)           // temp vector for matrix elements
    l = L   = J(n, 1, 1)           // indices of left boundary (old and new)
    r = R   = J(n, 1, q)           // indices of right boundary (old and new)
    nl      = 0                    // number of cells below left boundary
    nr      = n * q                // number of cells within right boundary
    ccw     = quadrunningsum(wneg) // cumulative column weights
    W       = quadsum(wpos:*ccw[rows(ccw)]) // sum of weights in search matrix
    Wl = WL = 0                    // sum of weights below left boundary
    WR      = W                    // sum of weights within right boundary
    k       = W/2                  // target quantile
    while ((nr-nl)>n) {
        // get trial value
        m = 0
        for (i=1; i<=n; i++) {
            if (l[i]<=r[i]) {
                // high median within row
                xx[++m] = -__mm_mc_el(xpos, xneg, npos, nzero, i,
                           l[i]+trunc((r[i]-l[i]+1)/2))
            }
        }
        trial = _mm_hl_qhi(xx[|1 \ m|], .5)
        // move right border
        j = 0
        for (i=n; i>=1; i--) {
            if (j<q) {
                while (-__mm_mc_el(xpos, xneg, npos, nzero, i, j+1)<trial) {
                    j++
                    if (j==q) break
                }
            }
            R[i] = j
        }
        p = (R:>0)
        if (any(p)) WR = quadsum(select(wpos, p) :* ccw[select(R, p)])
        else        WR = 0
        if (WR>k) {
            swap(r, R)
            nr = sum(R)
            continue
        }
        // move left border
        j = q + 1
        for (i=1; i<=n; i++) {
            while (-__mm_mc_el(xpos, xneg, npos, nzero, i, j-1)>trial) {
                j--
            }
            L[i] = j
        }
        p = (L:>1)
        WL = quadsum(select(wpos, p) :* ccw[select(L, p):-1])
        if (WL<k) {
            swap(l, L)
            Wl = WL
            nl = sum(L) - n
            continue
        }
        // trial = low quantile = high quantile
        if (WR==WL | (WR<k & WL>k)) return(-trial)
        // trial = low quantile
        if (WL==k) {
            m = 0
            for (i=1; i<=n; i++) {
                if (L[i]>q) continue
                xx[++m] = -__mm_mc_el(xpos, xneg, npos, nzero, i, L[i])
            }
            return(-(trial+min(xx[|1 \ m|]))/2)
        }
        // trial = high quantile
        for (i=1; i<=n; i++) {
            xx[i] = -__mm_mc_el(xpos, xneg, npos, nzero, i, R[i])
        }
        return(-(trial+max(xx))/2)
    }
    // get target value from remaining candidates
    m = 0
    for (i=1; i<=n; i++) {
        if (l[i]<=r[i]) {
            for (j=l[i]; j<=r[i]; j++) {
                m++
                xx[m] = -__mm_mc_el(xpos, xneg, npos, nzero, i, j)
                ww[m] = wpos[i] * wneg[j]
            }
        }
    }
    return(-_mm_hl_q_w(xx[|1 \ m|], ww[|1 \ m|], k, Wl))
}

real scalar _mm_mc_fw(real colvector x, real colvector w) // assumes med(x)=0
{
    real scalar     i, j, m, k, n, q, nl, nr, trial, npos, nzero, Wl, WR, WL, W
    real colvector  xx, ww, l, r, L, R, p, ccw, ccwz, xpos, xneg, wpos, wneg, wnegz, wzero

    if (any(trunc(w):!=w)) _error(3498, "non-integer frequency not allowed")
    if (rows(x)<=1) return(0)       // returning zero if n=1
    p = (x:>0)
    xpos    = select(x, p)          // obervations > median
    wpos    = select(w, p)          // weights of obervations >= median
    npos    = rows(xpos)            // number of obs > median
    p = (x:<0)
    xneg    = select(x, p)          // observations < median
    wneg    = select(w, p)          // weights of obervations <= median
    wzero   = sum(select(w, x:==0)) // aggregate weights for x==median
    nzero   = (wzero>0)             // has x==median
    if (nzero>0) {
        wpos  = wpos \ 1
        wnegz = (wzero>1 ? (comb(wzero, 2), wzero, comb(wzero, 2))' : wzero) \ wneg*wzero
        wneg  = wzero \ wneg
    }
    n       = npos + nzero          // number of rows in search matrix
    q       = nzero + rows(xneg)    // number of columns in search matrix
    xx = ww = J(n, 1, .)            // temp vector for matrix elements
    l = L   = J(n, 1, 1)            // indices of left boundary (old and new)
    r = R   = J(n, 1, q)            // indices of right boundary (old and new)
    if (wzero>1) {
        r[n] = q+2
        R[n] = q+2
    }
    nl      = 0                     // number of cells below left boundary
    nr      = sum(R)                // number of cells within right boundary
    ccw     = runningsum(wneg)      // cumulative column weights
    if (nzero==0) {
        W   = quadsum(wpos:*ccw[rows(ccw)]) // sum of weights in search matrix
    }
    else {
        ccwz = runningsum(wnegz)    // cumulative column weights in last row
        if (npos>=1) W = quadsum(wpos[|1 \ n-1|]:*ccw[rows(ccw)]) + ccwz[rows(ccwz)]
        else         W = ccwz[rows(ccwz)]
    }
    Wl = WL = 0                       // sum of weights below left boundary
    WR      = W                       // sum of weights within right boundary
    k       = W/2                     // target quantile
    while ((nr-nl)>n) {
        // get trial value
        m = 0
        for (i=1; i<=(n-(wzero>1)); i++) {
            if (l[i]<=r[i]) {
                // high median within row
                xx[++m] = -__mm_mc_el(xpos, xneg, npos, nzero, i,
                           l[i]+trunc((r[i]-l[i]+1)/2))
            }
        }
        if (wzero>1) { // handle last row
            if (l[n]<=r[n]) {
                m++
                xx[m] = -__mm_mc_el(xpos, xneg, npos-1, 3, n,
                           l[n]+trunc((r[n]-l[n]+1)/2))
            }
        }
        trial = _mm_hl_qhi(xx[|1 \ m|], .5)
        // move right border
        if (wzero>1) { // handle last row
            j = 0
            while (-__mm_mc_el(xpos, xneg, npos-1, 3, n, j+1)<trial) {
                j++
                if (j==(q+2)) break
            }
            R[n] = j
        }
        j = 0
        for (i=(n-(wzero>1)); i>=1; i--) {
            if (j<q) {
                while (-__mm_mc_el(xpos, xneg, npos, nzero, i, j+1)<trial) {
                    j++
                    if (j==q) break
                }
            }
            R[i] = j
        }
        p = (R:>0)
        if (nzero>0) p[n] = 0
        if (any(p)) WR = quadsum(select(wpos, p) :* ccw[select(R, p)])
        else        WR = 0
        if (nzero>0) {
            if (R[n]>0) WR = WR + ccwz[R[n]]
        }
        if (WR>k) {
            swap(r, R)
            nr = sum(R)
            continue
        }
        // move left border
        j = q + 1
        for (i=1; i<=(n-(wzero>1)); i++) {
            while (-__mm_mc_el(xpos, xneg, npos, nzero, i, j-1)>trial) {
                j--
            }
            L[i] = j
        }
        if (wzero>1) { // handle last row
            j = q + 3
            while (-__mm_mc_el(xpos, xneg, npos-1, 3, n, j-1)>trial) {
                j--
            }
            L[n] = j
        }
        p = (L:>1)
        if (nzero>0) p[n] = 0
        if (any(p)) WL = quadsum(select(wpos, p) :* 
            (rows(ccw)==1 ? ccw[select(L, p):-1]' : ccw[select(L, p):-1]))
        else        WL = 0
        if (nzero>0) {
            if (L[n]>1) WL = WL + ccwz[L[n]-1]
        }
        if (WL<k) {
            swap(l, L)
            Wl = WL
            nl = sum(L) - n
            continue
        }
        // trial = low quantile = high quantile
        if (ceil(k)!=k | (WR<k & WL>k)) return(-trial)
        // trial = low quantile
        if (WL==k) {
            m = 0
            for (i=1; i<=(n-(wzero>1)); i++) {
                if (L[i]>q) continue
                xx[++m] = -__mm_mc_el(xpos, xneg, npos, nzero, i, L[i])
            }
            if (wzero>1) { // handle last row
                if (L[n]<=(q+2)) {
                    xx[++m] = -__mm_mc_el(xpos, xneg, npos-1, 3, n, L[n])
                }
            }
            return(-(trial+min(xx[|1 \ m|]))/2)
        }
        // trial = high quantile
        for (i=1; i<=(n-(wzero>1)); i++) {
            xx[i] = -__mm_mc_el(xpos, xneg, npos, nzero, i, R[i])
        }
        if (wzero>1) { // handle last row
            xx[n] = -__mm_mc_el(xpos, xneg, npos-1, 3, n, R[n])
        }
        return(-(trial+max(xx))/2)
    }
    // get target value from remaining candidates
    m = 0
    for (i=1; i<=(n-nzero); i++) {
        if (l[i]<=r[i]) {
            for (j=l[i]; j<=r[i]; j++) {
                m++
                xx[m] = -__mm_mc_el(xpos, xneg, npos, nzero, i, j)
                ww[m] = wpos[i] * wneg[j]
            }
        }
    }
    if (nzero>0) { // handle last row
        if (l[n]<=r[n]) {
            for (j=l[n]; j<=r[n]; j++) {
                m++
                xx[m] = -__mm_mc_el(xpos, xneg, npos-(wzero>1),
                    (wzero>1 ? 3 : 1), n, j)
                ww[m] = wnegz[j]
            }
        }
    }
    return(-_mm_hl_q_w(xx[|1 \ m|], ww[|1 \ m|], k, Wl))
}

real scalar _mm_mc_naive(real colvector x) // noweights; assumes med(x)=0
{
    real scalar    i, j, m, n, q, npos, nzero
    real colvector xx, xpos, xneg 

    if (rows(x)<=1) return(0) // returning zero if n=1
    xpos   = select(x, x:>0)
    npos   = rows(xpos)
    xneg   = select(x, x:<0)
    nzero  = sum(x:==0)
    n      = npos + nzero
    q      = rows(xneg) + nzero
    m = 0
    xx = J(n*q, 1, .)
    for (i=1; i<=n; i++) {
        for (j=1; j<=q; j++) {
            xx[++m] = __mm_mc_el(xpos, xneg, npos, nzero, i, j)
        }
    }
    return(_mm_hl_q(xx, .5))
}

real scalar _mm_mc_naive_w(real colvector x, real colvector w) // assumes med(x)=0
{
    real scalar    i, j, m, n, q, npos, nzero
    real colvector xx, ww, xpos, xneg, wpos, wneg, wzero

    if (rows(x)<=1) return(0) // returning zero if n=1
    xpos   = select(x, x:>0)
    wpos   = select(w, x:>0)
    npos   = rows(xpos)
    xneg   = select(x, x:<0)
    wneg   = select(w, x:<0)
    wzero  = select(w, x:==0)
    nzero  = rows(wzero)
    if (nzero>0) {
        wpos = wpos \ wzero
        wneg = wzero[nzero::1] \ wneg // need to use reverse ordered wzero
    }
    n      = npos + nzero
    q      = rows(xneg) + nzero
    m = 0
    xx = ww = J(n*q, 1, .)
    for (i=1; i<=n; i++) {
        for (j=1; j<=q; j++) {
            m++
            xx[m] = __mm_mc_el(xpos, xneg, npos, nzero, i, j)
            ww[m] = wpos[i] * wneg[j]
        }
    }
    return(_mm_hl_q_w(xx, ww, .5))
}

real scalar _mm_mc_naive_fw(real colvector x, real colvector w) // assumes med(x)=0
{
    real scalar    i, j, m, n, q, npos, nzero, wzero
    real colvector xx, ww, xpos, xneg, wpos, wneg

    if (any(trunc(w):!=w)) _error(3498, "non-integer frequency not allowed")
    if (rows(x)<=1) return(0) // returning zero if n=1
    xpos   = select(x, x:>0)
    wpos   = select(w, x:>0)
    npos   = rows(xpos)
    xneg   = select(x, x:<0)
    wneg   = select(w, x:<0)
    wzero  = sum(select(w, x:==0)) // aggregate weights for x==median
    nzero  = (wzero>0)
    if (nzero>0) {
        wpos = wpos \ wzero
        wneg = wzero \ wneg
    }
    n      = npos + nzero
    q      = rows(xneg) + nzero
    m = 0
    xx = ww = J(n*q + 2*(wzero>1), 1, .)
    for (i=1; i<=n; i++) {
        for (j=1; j<=q; j++) {
            m++
            if (i>npos & j==nzero) { // x==median
                xx[m] = 0
                ww[m] = wzero
                if (wzero>1) {
                    m++
                    xx[m] = 1
                    ww[m] = comb(wzero, 2)
                    m++
                    xx[m] = -1
                    ww[m] = comb(wzero, 2)
                }
            }
            else {
                xx[m] = __mm_mc_el(xpos, xneg, npos, nzero, i, j)
                ww[m] = wpos[i] * wneg[j]
            }
        }
    }
    return(_mm_hl_q_w(xx, ww, .5))
}

// helper functions for quantiles

real scalar _mm_hl_q(real colvector x, real scalar P, 
    | real scalar offset) // must be integer; changes meaning of P if specified
{   // quantile (definition 2)
    real scalar    j0, j1, n, k
    real colvector p

    n = rows(x)
    if (n<1) return(.)
    if (n==1) return(x)
    if (args()==3) k = P     // P is a count (possibly noninteger)
    else           k = P * n // P is a proportion
    j0 = ceil(k)      - (args()==3 ? offset : 0) // index of low quantile
    j1 = floor(k) + 1 - (args()==3 ? offset : 0) // index of high quantile
    if (j0<1)      j0 = 1
    else if (j0>n) j0 = n
    if (j1<1)      j1 = 1
    else if (j1>n) j1 = n
    p = order(x, 1)
    if (j0==j1) return(x[p[j1]])
    return((x[p[j0]] + x[p[j1]])/2)
}

real scalar _mm_hl_qhi(real colvector x, real scalar P)
{   // high quantile
    real scalar    j, n
    real colvector p

    n = rows(x)
    if (n<1) return(.)
    p = order(x, 1)
    j = floor(P * n) + 1
    if (j<1)      j = 1
    else if (j>n) j = n
    return(x[p[j]])
}

real scalar _mm_hl_q_w(real colvector x, real colvector w, real scalar P, 
    | real scalar offset) // changes meaning of P if specified
{
    real scalar    n, i, k
    real colvector p, cw

    n = rows(x)
    if (n<1) return(.)
    p = order(x, 1)
    if (anyof(w, 0)) {
         p = select(p, w[p]:!=0)
         n = rows(p)
    }
    if (n<1) return(.)
    if (n==1) return(x[p])
    if (args()==4) {
        cw = quadrunningsum(offset \ w[p])[|2 \ n+1|]
        k = P // P is a count
    }
    else {
        cw = quadrunningsum(w[p])
        k = P * cw[n] // P is a proportion
    }
    if (k>=cw[n]) return(x[p[n]])
    for (i=1; i<=n; i++) {
        if (k>cw[i]) continue
        if (k==cw[i]) return((x[p[i]]+x[p[i+1]])/2)
        return(x[p[i]])
    }
    // cannot be reached
}

/*
real scalar _mm_hl_qhi_w(real colvector x, real colvector w, real scalar P)
{   // high quantile (weighted)
    real scalar    i, n, k
    real colvector p, cw

    n = rows(x)
    if (n<1) return(.)
    p = order(x, 1)
    cw = quadrunningsum(w[p])
    k  = cw[n] * P
    if (k>=cw[n]) return(x[p[n]])
    for (i=1; i<=n; i++) {
        if (k>=cw[i]) continue
        return(x[p[i]])
    }
    // cannot be reached
}
*/

end


*! {smcl}
*! {marker mm_ebalance}{bf:mm_ebalance.mata}{asis}
*! version 1.0.8  27apr2022  Ben Jann

version 11.2

// class & struct
local MAIN   mm_ebalance
local SETUP  _`MAIN'_setup
local Setup  struct `SETUP' scalar
local IF     _`MAIN'_IF
local If     struct `IF' scalar
// real
local RS     real scalar
local RR     real rowvector
local RC     real colvector
local RM     real matrix
// counters
local Int    real scalar
local IntC   real colvector
local IntR   real rowvector
local IntM   real matrix
// string
local SS     string scalar
// boolean
local Bool   real scalar
local BoolC  real colvector
// transmorphic
local T      transmorphic
local TS     transmorphic scalar
// pointers
local PC     pointer(real colvector) scalar
local PM     pointer(real matrix) scalar

mata:

// class ----------------------------------------------------------------------

struct `SETUP' {
    // data
    `PM'    X, X0      // pointers to main data and reference data
    `PC'    w, w0      // pointers to base weights
    `Int'   N, N0      // number of obs
    `RS'    W, W0      // sum of weights
    `RR'    m, m0      // means
    `RR'    s, s0      // scales
    `Int'   k          // number of terms
    `BoolC' omit       // flag collinear terms
    `Int'   k_omit     // number of omitted terms
    
    // settings
    `IntR'  adj, noadj // indices of columns to be adjusted/not adjusted
    `T'     tau        // target sum of weights
    `SS'    scale      // type of scales
    `RC'    btol       // balancing tolerance
    `SS'    ltype      // type of loss function
    `SS'    etype      // evaluator type
    `SS'    trace      // trace level
    `Int'   maxiter    // max number of iterations
    `RS'    ptol       // convergence tolerance for the parameter vector
    `RS'    vtol       // convergence tolerance for the balancing loss
    `Bool'  difficult  // use hybrid optimization
    `Bool'  nostd      // do not standardize
    `Bool'  nowarn     // do not display no convergence/balance warning
}

struct `IF' {
    `RM'   b, b0      // influence functions of coefficients
    `RC'   a, a0      // influence function of intercept
}

class `MAIN' {
    // settings
    private:
        void    new()             // initialize class with default settings
        void    clear()           // clear all results
        `Bool'  nodata()          // whether data is set
        `Setup' setup             // container for data and settings
    public:
        void    data()
        `RM'    X(), Xref()       // retrieve data
        `RC'    w(), wref()       // retrieve base weights
        `RS'    N(), Nref()       // retrieve number of obs
        `RS'    W(), Wref()       // retrieve sum of weights
        `RR'    m(), mref()       // retrieve moments
        `RR'    s(), sref()       // retrieve scales
        `RR'    mu()              // retrieve target moments
        `Int'   k()               // retrieve number of terms
        `RC'    omit()            // retrieve omitted flags
        `Int'   k_omit()          // retrieve number of omitted terms
        `T'     adj(), noadj()    // set/retrieve adj/noadj
        `T'     tau()             // set/retrieve target sum of weights
        `T'     scale()           // set/retrieve scales
        `T'     btol()            // set/retrieve balancing tolerance
        `T'     ltype()           // set/retrieve loss function
        `T'     etype()           // set/retrieve evaluator type
        `T'     alteval()         // set/retrieve alteval flag (old)
        `T'     trace()           // set/retrieve trace level
        `T'     maxiter()         // set/retrieve max iterations
        `T'     ptol()            // set/retrieve p-tolerance
        `T'     vtol()            // set/retrieve v-tolerance
        `T'     difficult()       // set/retrieve difficult flag
        `T'     nostd()           // set/retrieve nostd flag
        `T'     nowarn()          // set/retrieve nowarn flag
    
    // results
    public:
        `RC'    b()               // retrieve coefficients
        `RS'    a()               // retrieve normalizing intercept
        `RC'    xb()              // retrieve linear prediction
        `RC'    wbal()            // retrieve balancing weights
        `RC'    pr()              // retrieve propensity score
        `RR'    madj()            // adjusted (reweighted) means
        `RS'    wsum()            // retrieve sum of balancing weights
        `RS'    loss()            // retrieve final balancing loss
        `Bool'  balanced()        // retrieve balancing flag
        `RS'    value()           // retrieve value of optimization criterion
        `Int'   iter()            // retrieve number of iterations
        `Bool'  converged()       // retrieve convergence flag
        `RM'    IF_b(), IFref_b() // retrieve IF of coefficients
        `RC'    IF_a(), IFref_a() // retrieve IF of intercept
    private:
        `RS'    tau               // target sum of weight
        `RR'    mu                // target means
        `IntM'  adj, noadj        // permutation vectors for source of mu
        `RR'    scale             // scales for standardization
        `RC'    b                 // coefficients
        `RS'    a                 // normalizing intercept
        `RC'    xb                // linear prediction (without a)
        `RC'    wbal              // balancing weights
        `RR'    madj              // adjusted (reweighted) means
        `RS'    wsum              // sum of balancing weights
        `RC'    loss              // balancing loss
        `Bool'  balanced          // balance achieved
        `RS'    value             // value of optimization criterion
        `Int'   iter              // number of iterations
        `Bool'  conv              // optimize() convergence
        `If'    IF                // influence functions
        void    _IF_b(), _IF_a()  // generate influence functions
        void    Fit()             // fit coefficients
        void    _Fit_b(), _Fit_a()
        `RM'    _Fit_b_X(), _Fit_b_Xc()
        `RR'    _Fit_b_mu()
        void    _setadj()         // fill in adj and noadj
}

// init -----------------------------------------------------------------------

void `MAIN'::new()
{
    setup.tau       = "Wref"
    setup.scale     = "main"
    setup.ltype     = "reldif"
    setup.etype     = "bl"
    setup.nostd     = 0
    setup.trace     = (st_global("c(iterlog)")=="off" ? "none" : "value")
    setup.difficult = 0
    setup.maxiter   = st_numscalar("c(maxiter)")
    setup.ptol      = 1e-6
    setup.vtol      = 1e-7
    setup.btol      = 1e-6
    setup.nowarn    = 0
    setup.adj       = .
    b               = .z
}

void `MAIN'::clear()
{
    mu   = J(1,0,.)
    adj  = noadj = J(0,0,.)
    tau  = wsum = .
    b    = .z
    a    = .
    xb   = wbal = J(0,1,.)
    madj = J(1,0,.)
    loss = balanced = value = iter = conv = .
    IF   = `IF'()
}

// data -----------------------------------------------------------------------

void `MAIN'::data(`RM' X, `RC' w, `RM' X0, `RC' w0, | `Bool' fast)
{
    `RM' CP
    
    // check for missing values and negative weights
    if (args()<5) fast = 0
    if (!fast) {
        if (missing(X) | missing(w) | missing(X0) | missing(w0)) _error(3351)
        if (any(w:<0) | any(w0:<0)) _error(3498, "w and wref must be positive")
    }
    // obtain main data
    setup.k = cols(X)
    setup.X = &X
    setup.w = &w
    setup.N = rows(X)
    if (setup.N==0) _error(2000, "no observations in main data")
    if (rows(w)!=1) {
        if (rows(X)!=rows(w)) _error(3200, "X and w not conformable")
        setup.W = quadsum(w)
    }
    else setup.W = setup.N * w
    if (setup.N!=0 & setup.W==0) _error(3498, "sum(w) must be > 0")
    // obtain reference data
    setup.X0 = &X0
    setup.w0 = &w0
    setup.N0 = rows(X0)
    if (setup.N0==0) _error(2000, "no d observations in reference data")
    if (rows(w0)!=1) {
        if (rows(X0)!=rows(w0)) _error(3200, "Xref and wref not conformable")
        setup.W0 = quadsum(w0)
    }
    else setup.W0 = setup.N0 * w0
    if (setup.N0!=0 & setup.W0==0) _error(3498, "sum(wref) must be > 0")
    // if scale is set by user
    if (setup.scale=="user") {
        if (setup.k!=length(scale)) _error(3200, "X not conformable with scale")
    }
    else scale = J(1,0,.) // (clear scale)
    // target moments
    if (setup.k!=cols(X0)) _error(3200, "X and Xref not conformable")
    setup.m0 = mean(X0, w0)
    // identify collinear terms
    setup.m  = mean(X, w)
    if (setup.k==0) {
        setup.omit = J(0,1,.)
        setup.s = J(1,0,.)
        setup.k_omit = 0
    }
    else {
        CP = quadcrossdev(X, setup.m, w, X, setup.m)
        setup.s = sqrt(diagonal(CP)' / setup.W)
        setup.omit = (diagonal(invsym(CP)):==0) // or: diagonal(invsym(CP, 1..setup.k)):==0
        setup.k_omit = sum(setup.omit)
    }
    // clear results
    setup.s0 = J(1,0,.) // (scale of refdata will be set later only if needed)
    clear()
}

`Bool' `MAIN'::nodata() return(setup.X==NULL)

`RM' `MAIN'::X()
{
    if (nodata()) return(J(0,1,.))
    return(*setup.X)
}

`RM' `MAIN'::Xref()
{
    if (nodata()) return(J(0,1,.))
    return(*setup.X0)
}

`RC' `MAIN'::w()
{
    if (nodata()) return(J(0,1,.))
    return(*setup.w)
}

`RC' `MAIN'::wref()
{
    if (nodata()) return(J(0,1,.))
    return(*setup.w0)
}

`Int' `MAIN'::N() return(setup.N)

`Int' `MAIN'::Nref() return(setup.N0)

`RS' `MAIN'::W() return(setup.W)

`RS' `MAIN'::Wref() return(setup.W0)

`RR' `MAIN'::m() return(setup.m)

`RR' `MAIN'::mref() return(setup.m0)

`RR' `MAIN'::s() return(setup.s)

`RR' `MAIN'::sref() {
    if (length(setup.s0)) return(setup.s0)
    if (nodata()) return(setup.s0)
    setup.s0 = sqrt(diagonal(quadcrossdev(*setup.X0, setup.m0, *setup.w0, 
        *setup.X0, setup.m0))'/setup.W0)
    return(setup.s0)
}

`Int' `MAIN'::k() return(setup.k)

`RC'  `MAIN'::omit() return(setup.omit)

`Int' `MAIN'::k_omit() return(setup.k_omit)

// settings -------------------------------------------------------------------

`T' `MAIN'::tau(| `TS' tau)
{
    if (args()==0) {
        if (tau<.) return(tau)
        if (nodata()) return(setup.tau)
        if      (setup.tau=="Wref") tau = setup.W0
        else if (setup.tau=="W")    tau = setup.W
        else if (setup.tau=="Nref") tau = setup.N0
        else if (setup.tau=="N")    tau = setup.N
        else                        tau = setup.tau
        return(tau)
    }
    if (setup.tau==tau) return // no change
    if (isstring(tau)) {
        if (!anyof(("Wref", "W", "Nref", "N"), tau)) {
            printf("{err}'%s' not allowed\n", tau)
            _error(3498)
        }
    }
    else if (tau<=0 | tau>=.) _error(3498, "setting out of range")
    setup.tau = tau
    clear()
}

`T' `MAIN'::scale(| `T' scale0)
{
    `RR' s
    
    if (args()==0) {
        if (length(scale)) return(scale)
        if (setup.scale=="user") return(scale)
        if (nodata()) return(setup.scale)
        if      (setup.scale=="main") scale = s()
        else if (setup.scale=="ref")  scale = sref()
        else if (setup.scale=="avg")  scale = (s() + sref()) / 2
        else if (setup.scale=="wavg") scale = (s()*W() + sref()*Wref()) / 
                                              (W() + Wref())
        else if (setup.scale=="pooled") {
            scale = sqrt(diagonal(mm_variance0(X() \ Xref(), 
                (rows(w())==1  ? J(N(), 1, w())   : w()) \ 
                (rows(wref())==1 ? J(Nref(), 1, wref()) : wref())))')
        }
        _editvalue(scale, 0, 1)
        return(scale)
    }
    if (isstring(scale0)) {
        if (length(scale0)!=1) _error(3200)
        if (setup.scale==scale0) return // no change
        if (!anyof(("main", "ref", "avg", "wavg","pooled"), scale0)) {
            printf("{err}'%s' not allowed\n", scale0)
            _error(3498)
        }
        setup.scale = scale0
        scale = J(1,0,.)
    }
    else {
        s = editvalue(vec(scale0)', 0, 1)
        if (scale==s) return // no change
        if (missing(s)) _error(3351)
        if (any(s:<=0)) _error(3498, "scale must be positive")
        if (nodata()==0) {
            if (setup.k!=length(s)) _error(3200, "scale not conformable with X")
        }
        setup.scale = "user"
        scale = s
    }
    clear()
}

`T' `MAIN'::btol(| `RS' btol)
{
    if (args()==0) return(setup.btol)
    if (setup.btol==btol) return // no change
    if (btol<=0) _error(3498, "setting out of range")
    setup.btol = btol
    clear()
}

`T' `MAIN'::ltype(| `SS' ltype)
{
    if (args()==0) return(setup.ltype)
    if (setup.ltype==ltype) return // no change
    if (!anyof(("reldif", "absdif", "norm"), ltype)) {
        printf("{err}'%s' not allowed\n", ltype)
        _error(3498)
    }
    setup.ltype = ltype
    clear()
}

`T' `MAIN'::etype(| `SS' etype)
{
    if (args()==0) return(setup.etype)
    if (setup.etype==etype) return // no change
    if (!anyof(("bl","wl","mm","mma"), etype)) {
        printf("{err}'%g' not allowed\n", etype)
        _error(3498)
    }
    setup.etype = etype
    clear()
}

`T' `MAIN'::alteval(| `Bool' alteval) // for backward compatibility
{
    if (args()==0) return(setup.etype=="wl")
    if (setup.etype==(alteval!=0 ? "wl" : "bl")) return // no change
    setup.etype = (alteval!=0 ? "wl" : "bl")
    clear()
}

`T' `MAIN'::trace(| `SS' trace)
{
    `T' S
    
    if (args()==0) return(setup.trace)
    if (setup.trace==trace) return // no change
    S = optimize_init()
    optimize_init_tracelevel(S, trace) // throw error if trace is invalid
    setup.trace = trace
    clear()
}

`T' `MAIN'::maxiter(| `Int' maxiter)
{
    if (args()==0) return(setup.maxiter)
    if (setup.maxiter==maxiter) return // no change
    if (maxiter<0) _error(3498, "setting out of range")
    setup.maxiter = maxiter
    clear()
}

`T' `MAIN'::ptol(| `RS' ptol)
{
    if (args()==0) return(setup.ptol)
    if (setup.ptol==ptol) return // no change
    if (ptol<=0) _error(3498, "setting out of range")
    setup.ptol = ptol
    clear()
}

`T' `MAIN'::vtol(| `RS' vtol)
{
    if (args()==0) return(setup.vtol)
    if (setup.vtol==vtol) return // no change
    if (vtol<=0) _error(3498, "setting out of range")
    setup.vtol = vtol
    clear()
}

`T' `MAIN'::difficult(| `Bool' difficult)
{
    if (args()==0) return(setup.difficult)
    if (setup.difficult==(difficult!=0)) return // no change
    setup.difficult = (difficult!=0)
    clear()
}

`T' `MAIN'::nostd(| `Bool' nostd)
{
    if (args()==0) return(setup.nostd)
    if (setup.nostd==(nostd!=0)) return // no change
    setup.nostd = (nostd!=0)
    clear()
}

`T' `MAIN'::nowarn(| `Bool' nowarn)
{
    if (args()==0) return(setup.nowarn)
    if (setup.nowarn==(nowarn!=0)) return // no change
    setup.nowarn = (nowarn!=0)
    clear()
}

// target means ---------------------------------------------------------------

`T' `MAIN'::adj(| `RM' adj0)
{
    `RR' adj
    
    if (args()==0) {
        if (nodata()) return(setup.adj)
        _setadj()
        return(this.adj)
    }
    adj = mm_unique(trunc(vec(adj0)))'
    if (length(setup.noadj)==0) {
        if (setup.adj==adj) return // no change
    }
    if (adj!=.) {
        if (any(adj:<1)) _error(3498, "setting out of range")
    }
    setup.adj   = adj
    setup.noadj = J(1,0,.)
    clear()
}

`T' `MAIN'::noadj(| `RM' noadj0)
{
    `RR' noadj
    
    if (args()==0) {
        if (nodata()) return(setup.noadj)
        _setadj()
        return(this.noadj)
    }
    noadj = mm_unique(trunc(vec(noadj0)))'
    if (setup.adj==.) {
        if (setup.noadj==noadj) return // no change
    }
    if (any(noadj:<1)) _error(3498, "setting out of range")
    setup.adj   = .
    setup.noadj = noadj
    clear()
}

void `MAIN'::_setadj() // assumes that data has been set
{
    `Int'  i, j, k
    `IntR' p
    
    if (rows(adj)) return // already set
    // zero variables
    k = k()
    if (k==0) {
        adj = noadj = J(1,0,.)
        return
    }
    // case 1: noadj() has been set
    if (length(setup.noadj)) {
        p = J(1,k,1)
        for (i=length(setup.noadj); i; i--) {
            j = setup.noadj[i]
            if (j>k) continue // be tolerant and ignore invalid subscripts
            p[j] = 0
        }
    }
    // case 2: adj() has been set
    else if (setup.adj!=.) {
        p = J(1,k,0)
        for (i=length(setup.adj); i; i--) {
            j = setup.adj[i]
            if (j>k) continue // be tolerant and ignore invalid subscripts
            p[j] = 1
        }
    }
    // case 3: default (adjust all)
    else p = 1
    // fill in adj/noadj
    if (allof(p, 1)) {
        adj   = .
        noadj = J(1,0,.)
    }
    else {
        adj   = select(1..k, p)
        noadj = select(1..k,!p)
    }
}

`RR' `MAIN'::mu()
{
    if (length(mu)) return(mu)
    if (nodata())   return(mu)
    if (adj()==.) {
        mu = setup.m0
    }
    else {
        mu = setup.m
        if (length(adj())) mu[adj()] = setup.m0[adj()]
    }
    return(mu)
}

// results --------------------------------------------------------------------

`RC' `MAIN'::b()
{
    if (b==.z) Fit()
    return(b)
}

`RS' `MAIN'::a()
{
    if (b==.z) Fit()
    return(a)
}

`RC' `MAIN'::wbal()
{
    if (b==.z) Fit()
    return(wbal)
}

`RC' `MAIN'::xb()
{
    if (b==.z) Fit()
    return(xb :+ a)
}

`RC' `MAIN'::pr()
{
    if (b==.z) Fit()
    return(invlogit(xb :+ (a + ln(Wref()/tau()))))
}

`RR' `MAIN'::madj()
{
    if (length(madj)) return(madj)
    if (b==.z) Fit()
    madj = mean(X(), wbal)
    return(madj)
}

`RS' `MAIN'::wsum()
{
    if (wsum<.) return(wsum)
    if (b==.z) Fit()
    wsum = quadsum(wbal)
    return(wsum)
}

`RS' `MAIN'::loss()
{
    if (b==.z) Fit()
    return(loss)
}

`RS' `MAIN'::balanced()
{
    if (b==.z) Fit()
    return(balanced)
}

`RS' `MAIN'::value()
{
    if (b==.z) Fit()
    return(value)
}

`Int' `MAIN'::iter()
{
    if (b==.z) Fit()
    return(iter)
}

`Bool' `MAIN'::converged()
{
    if (b==.z) Fit()
    return(conv)
}

`RM' `MAIN'::IF_b()
{
    if (rows(IF.b)==0) _IF_b()
    return(IF.b)
}

`RM' `MAIN'::IFref_b()
{
    if (rows(IF.b0)==0) _IF_b()
    return(IF.b0)
}

`RC' `MAIN'::IF_a()
{
    if (rows(IF.a)==0) _IF_a()
    return(IF.a)
}

`RC' `MAIN'::IFref_a()
{
    if (rows(IF.a0)==0) _IF_a()
    return(IF.a0)
}

// optimization ---------------------------------------------------------------

void `MAIN'::Fit()
{
    `Bool' nofit
    
    // optimize
    if (nodata()) _error(3498, "data not set")
    nofit = 0
    if ((k()-k_omit())<=0)   nofit = 1     // no covariates
    else if (!length(adj())) nofit = 1     // no adjustments
    else if (length(noadj()) & k_omit()) { // no adjustments among non-omitted
        nofit = all(omit()[adj()])
    }
    if (nofit) {
        b = J(k(), 1, 0)
        iter = value = 0
        conv = 1
    }
    else _Fit_b() // fit coefficients
    _Fit_a()      // compute intercept (and balancing weights)
    
    // check balancing
    if (etype()!="bl" | nostd()==0 | k_omit()) {
        // compute balancing loss using raw data
        loss = _mm_ebalance_loss(ltype(), mean(X():-mu(), wbal), mu())
    }
    else loss = value
    if (trace()!="none") {
        printf("{txt}Final fit:     balancing loss = {res}%10.0g\n", loss)
    }
    balanced = (loss<btol())
    if (balanced) return
    if (nowarn()) return
    display("{err}balance not achieved")
}

void `MAIN'::_Fit_a()
{
    `RS' ul
    
    xb = X() * b
    ul = max(xb) // set exp(max)=1 to avoid numerical overflow
    a  = ln(tau()) - ln(quadsum(w() :* exp(xb :- ul))) - ul
    wbal = w() :* exp(xb :+ a)
}

void `MAIN'::_Fit_b()
{
    `RR'   beta
    `IntC' p
    `T'    S
    
    // setup
    if (k_omit()) p = select(1::k(), omit():==0)
    S = optimize_init()
    optimize_init_which(S, "min")
    optimize_init_technique(S, "nr")
    optimize_init_evaluatortype(S, "d2")
    optimize_init_conv_ignorenrtol(S, "on")
    optimize_init_tracelevel(S, trace())
    optimize_init_conv_maxiter(S, maxiter())
    optimize_init_conv_ptol(S, ptol())
    optimize_init_conv_vtol(S, vtol())
    optimize_init_singularHmethod(S, difficult() ? "hybrid" : "")
    optimize_init_conv_warning(S, nowarn() ? "off" : "on")
    if (etype()=="mma") {
        optimize_init_evaluator(S, &_mm_ebalance_mma())   // gmm w/ alpha
        optimize_init_valueid(S, "criterion Q(p)")
        optimize_init_params(S, J(1, k()-k_omit()+1, 0)) // starting values
        optimize_init_argument(S, 1, _Fit_b_X(p))        // data
        optimize_init_argument(S, 2, _Fit_b_mu(p))       // target moments
        optimize_init_argument(S, 3, w()/W())            // norm. base weights
        optimize_init_argument(S, 4, tau()/W())          // normalized tau
    }
    else if (etype()=="mm") {
        optimize_init_evaluator(S, &_mm_ebalance_mm())   // gmm
        optimize_init_valueid(S, "criterion Q(p)")
        optimize_init_params(S, J(1, k()-k_omit(), 0))   // starting values
        optimize_init_argument(S, 1, _Fit_b_X(p))        // data
        optimize_init_argument(S, 2, _Fit_b_mu(p))       // target moments
        optimize_init_argument(S, 3, w())                // base weights
    }
    else if (etype()=="wl") {
        optimize_init_evaluator(S, &_mm_ebalance_lw())   // sum(ln(w))
        optimize_init_valueid(S, "criterion L(w)")
        optimize_init_params(S, J(1, k()-k_omit(), 0))   // starting values
        optimize_init_argument(S, 1, _Fit_b_Xc(p))       // centered data
        optimize_init_argument(S, 2, w())                // base weights
    }
    else {
        optimize_init_evaluator(S, &_mm_ebalance_bl())  // balance loss
        optimize_init_valueid(S, "balancing loss")
        optimize_init_params(S, J(1, k()-k_omit(), 0))  // starting values
        optimize_init_argument(S, 1, _Fit_b_Xc(p))      // centered data
        optimize_init_argument(S, 2, _Fit_b_mu(p))      // target moments
        optimize_init_argument(S, 3, w())               // base weights
        optimize_init_argument(S, 4, ltype())           // loss type
    }
    
    // run optimizer
    (void) _optimize(S)
    if (optimize_result_errorcode(S)) {
        errprintf("{p}\n")
        errprintf("%s\n", optimize_result_errortext(S))
        errprintf("{p_end}\n")
        exit(optimize_result_returncode(S))
    }
    
    // obtain results
    beta = optimize_result_params(S)
    if (etype()=="mma") beta = beta[|1\length(beta)-1|] // discard alpha
    if (k_omit()) {
        b = J(k(), 1, 0)
        if (nostd()) b[p] = beta'
        else         b[p] = (beta :/ scale()[p])'
    }
    else {
        if (nostd()) b = beta'
        else         b = (beta :/ scale())'
    }
    iter  = optimize_result_iterations(S)
    value = optimize_result_value(S)
    conv  = optimize_result_converged(S)
}

`RM' `MAIN'::_Fit_b_X(`IntC' p)
{
    if (k_omit()) {
        if (nostd()) return(X()[,p])
        return(X()[,p] :/ scale()[p])
    }
    if (nostd()) return(X())
    return(X() :/ scale())
}

`RM' `MAIN'::_Fit_b_Xc(`IntC' p)
{
    if (k_omit()) {
        if (nostd()) return(X()[,p] :- mu()[p])
        return((X()[,p] :- mu()[p]) :/ scale()[p])
    }
    if (nostd()) return(X() :- mu())
    return((X() :- mu()) :/ scale())
}

`RR' `MAIN'::_Fit_b_mu(`IntC' p)
{
    if (k_omit()) {
        if (nostd()) return(mu()[p])
        return((mu() :/ scale())[p])
    }
    if (nostd()) return(mu())
    return(mu() :/ scale())
}

`RS' _mm_ebalance_loss(`SS' ltype, `RR' d, `RR' mu)
{
    if (ltype=="absdif") return(max(abs(d)))
    if (ltype=="norm") return(sqrt(d*d'))
    return(mreldif(d+mu, mu)) // ltype=="reldif"
}

void _mm_ebalance_bl(`Int' todo, `RR' b, `RM' X, `RR' mu, `RC' w0, `SS' ltype,
    `RS' v, `RR' g, `RM' H)
{   // evaluator based on balance loss
    `RS' W
    `RC' w
    
    w = X * b'
    w = w0 :* exp(w :- max(w)) // avoid numerical overflow
    W = quadsum(w)
    g = quadcross(w, X) / W
    v = _mm_ebalance_loss(ltype, g, mu)
    if (todo==2) H = quadcross(X, w, X) / W
}

void _mm_ebalance_lw(`Int' todo, `RR' b, `RM' X, `RC' w0,
    `RS' v, `RR' g, `RM' H)
{   // evaluator using sum(ln(weights)) as criterion
    `RS' W
    `RC' w
    
    w = X * b'
    W = max(w)
    w = w0 :* exp(w :- W) // avoid numerical overflow
    v = ln(quadsum(w)) + W
    if (todo>=1) {
        W = quadsum(w)
        g = quadcross(w, X) / W
        if (todo==2) H = quadcross(X, w, X) / W
    }
}

void _mm_ebalance_mm(`Int' todo, `RR' b, `RM' X, `RR' mu, `RC' w0,
    `RS' v, `RR' g, `RM' H)
{   // gmm type evaluator
    `RC' w
    `RR' d
    `RM' h, G
    
    w = X * b'
    w = w0 :* exp(w :- max(w))   // avoid numerical overflow
    w = w :/ quadcolsum(w)
    h = w :* (X :- mu)
    d = quadcolsum(h)
    v = d * d'
    if (todo>=1) {
        G = quadcross(h, X :- quadcolsum(w:*X))
        g = d * G
        if (todo==2) H = G'G 
    }
}

void _mm_ebalance_mma(`Int' todo, `RR' b, `RM' X, `RR' mu, `RC' w0, `RS' tau, 
    `RS' v, `RR' g, `RM' H)
{   // gmm type evaluator including alpha in optimization problem
    `RS' d_a
    `RC' w, h_a
    `RR' d_b
    `RM' h_b, G
    
    w   = w0 :* exp(X*b[|1\cols(b)-1|]' :+ b[cols(b)])
    h_b = w :* (X :- mu)
    h_a = w :- w0:*tau  // = w0 * (e(xb+a) - tau/sum(w0))
    d_b = quadcolsum(h_b)
    d_a = quadsum(h_a)
    v   = (d_b, d_a) * (d_b, d_a)'
    if (todo>=1) {
        G = (quadcross(h_b, X), d_b') \ (quadcross(w, X), colsum(w))
        g = (d_b, d_a) * G
        if (todo==2) H = G'G 
    }
}

// influence functions --------------------------------------------------------

void `MAIN'::_IF_b()
{
    if (b==.z) Fit()
    _mm_ebalance_IF_b(IF, X(), Xref(), w(), wbal, madj(), mu(), tau(), W(),
        Wref(), adj(), noadj(), omit())
}

void `MAIN'::_IF_a()
{
    if (rows(IF.b)==0) _IF_b()
    _mm_ebalance_IF_a(IF, X(), w(), wbal, tau(), W())
}

void _mm_ebalance_IF_b(`If' IF, `RM' X, `RM' Xref, `RC' w, `RC' wbal, `RR' madj,
    `RR' mu, `RS' tau, `RS' W, `RS' Wref, `IntR' adj, `IntR' noadj,
    | `BoolC' omit)
{   // using "alternative approach" formulas
    `Int'  k
    `IntC' p
    `RM'   G, Q
    
    G = quadcrossdev(X, mu, wbal/tau, X, madj)
    if (length(omit)==0) {
        // no information on omitted terms; use qrinv()
        Q = -qrinv(G)
    }
    else if (any(omit)) {
        // discard omitted terms during inversion
        k = cols(X)
        p = select(1::k, omit:==0)
        Q = J(k, k, 0)
        if (length(p)) Q[p,p] = -luinv(G[p,p])
    }
    else {
        // no omitted terms; save to use luinv()
        Q = -luinv(G)
    }
    IF.b  = (wbal/tau):/w :* (X :- madj) * Q' // [sic!] using madj instead of mu
    if (adj==.) IF.b0 = (-1/Wref) * (Xref :- mu) * Q'
    else {
        IF.b0 = J(rows(Xref), cols(Xref), 0)
        if (length(adj)) IF.b0[, adj] = 
            (-1/Wref) * (Xref[,adj] :- mu[adj]) * Q[adj,adj]'
        if (length(noadj)) IF.b[,noadj] = 
            IF.b[,noadj] - (X[,noadj] :- mu[noadj])/W * Q[noadj,noadj]'
    }
}

void _mm_ebalance_IF_a(`If' IF, `RM' X, `RC' w, `RC' wbal, `RS' tau, `RS' W)
{    // using simplified IF assuming tau as fixed
    `RM' Q
    
    Q = quadcross(wbal, X)
    IF.a  = ((wbal:/w :- tau/W) + IF.b * Q') / -tau
    IF.a0 = (IF.b0 * Q') / -tau
}

end



*! {smcl}
*! {marker mm_wbal}{bf:mm_wbal.mata}{asis}
*! version 1.0.0  25jul2021  Ben Jann
version 11.2

mata:

real colvector mm_wbal(real matrix X, real colvector w, 
    real matrix X0, real colvector w0, | real scalar nowarn)
{
    class mm_ebalance scalar S
    
    if (args()<5) nowarn = 0
    S.nowarn(nowarn)
    S.trace("none")
    S.data(X, w, X0, w0)
    return(S.wbal())
}

end



*! {smcl}
*! {marker mm_ebal}{bf:mm_ebal.mata}{asis}
*! version 1.0.4  05aug2019  Ben Jann
version 11.2
local Bool  real scalar
local BoolR real rowvector
local Int   real scalar
local IntV  real vector
local IntC  real colvector
local RS    real scalar
local RV    real vector
local RC    real colvector
local RR    real rowvector
local RM    real matrix
local SS    string scalar
local T     transmorphic
local pRC   pointer (`RC') scalar
local pRM   pointer (`RM') scalar
local S     struct _mm_ebal_struct scalar
mata:

struct _mm_ebal_struct {
    `T'    O         // optimization object
    `RS'   N         // target group size (sum of weights)
    `RM'   C         // constraint matrix (excluding collinear columns)
    `RM'   CC        // collinear columns from constraint matrix
    `RC'   Q         // base weights
    `RC'   W         // balancing weights
    `Bool' nc        // normalizing constraint included in optimization
    `RS'   btol      // balancing tolerance
    `Bool' balanced  // 1 if balance achieved, 0 else
    `RR'   Z         // coefficients
    `RR'   g         // gradient
    `RS'   v         // max difference
    `Bool' conv      // 1 if optimize() converged
    `Int'  rc        // return code of optimize()
    `Int'  i         // number of iterations
}

`RR'   mm_ebal_N(`S' S)        return(S.N)
`RM'   mm_ebal_C(`S' S)        return(S.C)
`RM'   mm_ebal_CC(`S' S)       return(S.CC)
`RC'   mm_ebal_Q(`S' S)        return(S.Q)
`RC'   mm_ebal_W(`S' S)        return(S.W)
`Bool' mm_ebal_balanced(`S' S) return(S.balanced)
`RR'   mm_ebal_g(`S' S)        return(S.g)
`RS'   mm_ebal_v(`S' S)        return(S.v)
`Bool' mm_ebal_conv(`S' S)     return(S.conv)
`Int'  mm_ebal_rc(`S' S)       return(S.rc)
`Int'  mm_ebal_i(`S' S)        return(S.i)

`S' mm_ebal_init(
    `RM'   X1,       // treatment group data
    `RC'   w1,       // treatment group base weights
    `RM'   X0,       // control group data
    `RC'   w0,       // control group base weights
  | `IntV' tar,      // balancing targets
    `Bool' cov,      // whether to balance covariances
    `Bool' nc,       // include normalizing constraint in optimization
    `Bool' dfc,      // apply degrees of freedom correction
    `Bool' nostd)    // do not standardize data
{
    `S'     S        // main struct
    `RM'    C1       // expanded X1
    `RM'    M        // target moments
    `RR'    sd       // standard deviations of C1
    `BoolR' collin   // collinear terms
    `RS'    N0       // size of control group (sum of weights)
    
    // defaults
    if (args()<5)  tar   = 1
    if (args()<6)  cov   = 0
    if (args()<7)  nc    = 0
    if (args()<8)  dfc   = 0
    if (args()<9)  nostd = 0
    S.nc   = nc
    S.btol = 1e-5
    // compute group sizes and check input
    // - treatment group
    if (hasmissing(X1)) _error(3351, "missing values in X1 not allowed")
    if (hasmissing(w1)) _error(3351, "missing values in w1 not allowed")
    if (any(w1:<0))     _error(3498, "negative values in w1 not allowed")
    if (rows(w1)==1) S.N = w1*rows(X1)
    else {
        if (rows(w1)!=rows(X1)) _error(3200, "X1 and w1 not conformable")
        S.N = quadsum(w1)
    }
    if (S.N<=0) _error(3498, "treatment group size (sum of weights) must be > 0")
    // - control group
    if (hasmissing(X0)) _error(3351, "missing values in X0 not allowed")
    if (hasmissing(w0)) _error(3351, "missing values in w0 not allowed")
    if (any(w0:<0))     _error(3498, "negative values in w0 not allowed")
    if (rows(w0)==1) N0 = w0*rows(X0)
    else {
        if (rows(w0)!=rows(X0)) _error(3200, "X0 and w0 not conformable")
        N0 = quadsum(w0)
    }
    if (N0<=0) _error(3498, "control group size (sum of weights) must be > 0")
    // - number of variables
    if (cols(X1)!=cols(X0)) _error(3200, "X1 and X0 must contain the same number of columns")
    // expand X1 and X0
    C1  = _mm_ebal_C(tar, cov, X1)
    S.C = _mm_ebal_C(tar, cov, X0)
    // determine collinear columns (across both groups)
    collin = _mm_ebal_collin((C1 \ S.C),
             ((rows(w1)==1 ? J(rows(C1), 1, w1) : w1)
            \ (rows(w0)==1 ? J(rows(S.C), 1, w0) : w0)))
    // compute target moments
    M = mean(C1, w1)
    // apply degrees-of-freedom correction to target moments
    if (dfc) _mm_ebal_dfc(tar, cov, M, cols(X1), (1-1/S.N) / (1-1/N0))
    // remove collinear columns
    if (any(collin)) {
        C1  = select(C1, collin:==0)
        M   = select(M, collin:==0)
        S.C = select(S.C, collin:==0)
    }
    // determine scale for standardization
    if (rows(C1)==1) sd = J(1, cols(C1), 1)
    else {
        if (nostd | cols(C1)==0) sd = J(1, cols(C1), 1)
        else {
            sd = sqrt(diagonal(variance(C1, w1)))'
            _editvalue(sd, 0, 1) // (sd must not be zero)
        }
    }
    // prepare constraints matrix
    // - determine collinear columns
    collin = _mm_ebal_collin(S.C, w0)
    // - standardize constraints
    S.C = (S.C :- M) :/ sd
    // - put collinear columns aside
    if (any(collin)) {
        S.CC = select(S.C, collin)
        S.C  = select(S.C, collin:==0)
    }
    // - add normalization constraint
    if (S.nc) S.C = J(rows(X0),1,1), S.C
    // prepare base weights
    if (S.nc) S.Q = w0
    else      S.Q = w0 / N0
    // prepare optimization object
    S.O = optimize_init()
    optimize_init_which(S.O, "min")
    optimize_init_evaluator(S.O, &_mm_ebal_eval())
    optimize_init_evaluatortype(S.O, "d2")
    optimize_init_conv_ignorenrtol(S.O, "on")
    optimize_init_valueid(S.O, "max difference")
    // set starting values
    if (S.nc) S.Z = (ln(S.N/N0), J(1, cols(S.C)-1, 0))
    else      S.Z = J(1, cols(S.C), 0)
    // done
    return(S)
}

`RM' _mm_ebal_C(`RV' t, `Bool' cov, `RM' X)
{
    `RS' i, j, k, l, c
    `RM' C
    
    // confirm that t only contains 1, 2, or 3
    l = length(t)
    for (i=1; i<=l; i++) {
        if (!anyof((1,2,3), t[i])) _error("invalid target")
    }
    // means
    k = cols(X)
    C = X
    // variances
    if (any(t:>=2)) {
        c = 0
        for (i=1; i<=k; i++) { // count number of terms
            if (t[mod(i-1, l) + 1] >= 2) c++
        }
        C = C, J(rows(C), c, .)
        c = cols(C) - c
        for (i=1; i<=k; i++) {
            if (t[mod(i-1, l) + 1] >= 2) {
                C[,++c] = X[,i]:^2
            }
        }
    }
    // covariances
    if (cov) {
        c = (k^2 - k)/2
        C = C, J(rows(C), c, .)
        c = cols(C) - c
        for (i=1; i<k; i++) {
            for (j=i+1; j<=k; j++) {
                C[,++c] = X[,i] :* X[,j]
            }
        }
    }
    // skewnesses
    if (anyof(t,3)) {
        c = 0
        for (i=1; i<=k; i++) { // count number of terms
            if (t[mod(i-1, l) + 1] == 3) c++
        }
        C = C, J(rows(C), c, .)
        c = cols(C) - c
        for (i=1; i<=k; i++) {
            if (t[mod(i-1, l) + 1] == 3) {
                C[,++c] = X[,i]:^3
            }
        }
    }
    return(C)
}

`BoolR' _mm_ebal_collin(`RM' C, `RC' w)
{
    `RR' m
    `RM' CP

    if (cols(C)==0) return(J(1,0,.))
    m  = mean(C, w)
    CP = quadcrossdev(C, m, w, C, m)
    return((diagonal(invsym(CP, 1..cols(CP))):==0)')
}

void _mm_ebal_dfc(`RV' t, `Bool' cov, `RM' M, `RS' k, `RS' dfc)
{
    `RS' i, j, l, c
    `RM' V
    
    if (cols(M)==k) return
    if (dfc==1)     return
    c = k
    l = length(t)
    // variances
    if (any(t:>=2)) {
        V = J(2, k, 1) // needed for skewness correction
        for (i=1; i<=k; i++) {
            if (t[mod(i-1, l) + 1] >= 2) {
                c++
                V[1,i] = M[c]
                M[c]   = (M[c] - (1-dfc) * M[i]^2) / dfc
                V[2,i] = M[c]
            }
        }
    }
    // covariances
    if (cov) {
        for (i=1; i<k; i++) {
            for (j=i+1; j<=k; j++) {
                c++
                M[c] = (M[c] - (1-dfc) * M[i] * M[j]) / dfc
            }
        }
    }
    // skewnesses
    if (anyof(t,3)) {
        for (i=1; i<=k; i++) {
            if (t[mod(i-1, l) + 1] == 3) {
                c++
                M[c] = (M[c] - (V[1,i] - dfc^1.5 * V[2,i]) * 3 *M[i]
                             + (1 - dfc^1.5) * 2 * M[i]^3) / dfc^1.5
            }
        }
    }
}

`T' mm_ebal_btol(`S' S, | `RS' btol)
{
    if (args()==1) return(S.btol)
    S.btol = btol
}

`T' mm_ebal_trace(`S' S, | `SS' tracelevel)
{
    if (args()==1) return(optimize_init_tracelevel(S.O))
    optimize_init_tracelevel(S.O, tracelevel)
}

`T' mm_ebal_difficult(`S' S, | `Bool' flag)
{
    if (args()==1) return(optimize_init_singularHmethod(S.O)=="hybrid")
    if (flag)      optimize_init_singularHmethod(S.O, "hybrid")
    else           optimize_init_singularHmethod(S.O, "")
}

`T' mm_ebal_maxiter(`S' S, | `Int' maxiter)
{
    if (args()==1) return(optimize_init_conv_maxiter(S.O))
    optimize_init_conv_maxiter(S.O, maxiter)
}

`T' mm_ebal_ptol(`S' S, | `RS' tol)
{
    if (args()==1) return(optimize_init_conv_ptol(S.O))
    optimize_init_conv_ptol(S.O, tol)
}

`T' mm_ebal_vtol(`S' S, | `RS' tol)
{
    if (args()==1) return(optimize_init_conv_vtol(S.O))
    optimize_init_conv_vtol(S.O, tol)
}

`T' mm_ebal_nowarn(`S' S, | `Bool' flag)
{
    if (args()==1) return(optimize_init_conv_warning(S.O)=="off")
    if (flag)      optimize_init_conv_warning(S.O, "off")
    else           optimize_init_conv_warning(S.O, "on")
}

`T' mm_ebal_Z(`S' S, | `RR' Z)
{
    if (args()==1) return(S.Z)
    S.Z = Z
}

`Bool' mm_ebal(`S' S)
{
    // find solution
    if (S.nc==0 & cols(S.C)==0) {
        // no covariates; nothing to do
        S.Z = S.g = J(1,0,.)
        S.conv = 1
        S.rc = S.i = S.v = 0
    }
    else {
        optimize_init_argument(S.O, 1, S.nc)
        optimize_init_argument(S.O, 2, &S.Q)
        optimize_init_argument(S.O, 3, &S.C)
        optimize_init_argument(S.O, 4, S.N)
        optimize_init_params(S.O, S.Z)
        (void) _optimize(S.O)
        if (optimize_result_errorcode(S.O)) {
            errprintf("{p}\n")
            errprintf("%s\n", optimize_result_errortext(S.O))
            errprintf("{p_end}\n")
        }
        S.Z    = optimize_result_params(S.O)
        S.conv = optimize_result_converged(S.O)
        S.rc   = optimize_result_returncode(S.O)
        S.i    = optimize_result_iterations(S.O)
        S.v    = optimize_result_value(S.O)
        S.g    = optimize_result_gradient(S.O)
    }
    // obtain balancing weights from solution and update balancing loss
    if (S.nc) {
        S.W = S.Q :* exp(quadcross(S.C', S.Z'))
        if (cols(S.CC)>0) S.v = max((S.v, abs(quadcross(S.W, S.CC)) / S.N))
    }
    else {
        S.W = quadcross(S.C', S.Z')
        S.W = S.Q :* exp(S.W:-max(S.W))
        S.W = S.W / quadsum(S.W)
        if (cols(S.CC)>0) S.v = max((S.v, abs(quadcross(S.W, S.CC))))
        S.W = S.W * S.N
    }
    if (cols(S.CC)>0) {
        if (optimize_init_tracelevel(S.O)!="none") {
            printf("{txt}Final fit including collinear terms:\n")
            printf("               max difference = {res}%10.0g\n", S.v)
        }
    }
    // check whether balancing criterion fulfilled
    S.balanced = (S.v < S.btol)
    return(S.balanced)
}

void _mm_ebal_eval(`Int' todo, `RR' Z, 
    `Bool' nc, `pRC' Q, `pRM' C, `RS' N, 
    `RS' v, `RR' g, `RM' H)
{
    `RC' W
    
    if (nc) {
        W = *Q :* exp(quadcross(*C', Z'))
        g = quadcross(W, *C)
        g[1] = g[1] - N
        v = max(abs(g)) / N
    }
    else {
        W = quadcross(*C', Z')
        W = *Q :* exp(W:-max(W)) // set exp(max)=1 to avoid numerical overflow
        W = W / quadsum(W)
        g = quadcross(W, *C)
        v = max(abs(g))
    }
    if (todo==2) {
        H = quadcross(*C, W, *C)
    }
}

end



*! {smcl}
*! {marker mm_crosswalk}{bf:mm_crosswalk.mata}{asis}
*! version 1.0.0  Ben Jann  23aug2021

version 11.2
mata:

transmorphic vector mm_crosswalk(
    transmorphic vector x,
    transmorphic vector from,
    transmorphic vector to,
    | transmorphic vector d, 
      transmorphic scalar n)
{
    real scalar         usehash, l, offset
    real rowvector      minmax
    transmorphic vector res

    // defaults
    if (args()<4) d = missingof(to)
    if (args()<5) n = 1e6
    
    // consistency checks
    if (eltype(x)!=eltype(from))  _error(3250, "{it:x} and {it:from} must be of same type")
    if (length(from)!=length(to)) _error(3200, "{it:from} and {it:to} must have same length")
    if (eltype(to)!=eltype(d))    _error(3250, "{it:to} and {it:d} must be of same type")
    if (length(d)!=1) {
        if (length(d)!=length(x)) _error(3200, "{it:x} and {it:d} not conformable")
    }
    
    // nothing to translate
    // - empty x
    if (!length(x)) return(J(rows(x), cols(x), missingof(to)))
    // - empty dictionary
    if (!length(to)) {
        if (length(d)!=1) return(J(1,1, rows(x)==rows(d) ? d : d'))
        return(J(rows(x), cols(x), d))
    }
    
    // pick algorithm
    usehash = 0
    if (!isreal(x)) usehash = 1
    else if (n<1)   usehash = 1
    else {
        if (n<.z) {
            if      (hasmissing(from))  usehash = 1
            else if (hasmissing(x))     usehash = 1
            else if (trunc(from)!=from) usehash = 1
            else if (trunc(x)!=x)       usehash = 1
        }
        if (!usehash) {
            minmax = minmax((minmax(from),minmax(x)))
            l = minmax[2] - minmax[1] + 1
            if (l>n) usehash = 1
        }
    }
    
    // use hash algorithm
    if (usehash) return(_mm_crosswalk_hash(x, from, to, d))
    
    // use index algorithm
    offset = minmax[1] - 1
    res = rows(to)!=1 ? J(l, 1, d[1]) : J(1, l, d[1])
    if (length(d)!=1) res[x :- offset] = (rows(res)!=1)==(rows(d)!=1) ? d : d'
    res[from :- offset] = to
    res = res[x :- offset]
    return(rows(res)==rows(x) ? res : res')
}

transmorphic vector mm_crosswalk_hash(
    transmorphic vector x,
    transmorphic vector from,
    transmorphic vector to,
    | transmorphic vector d)
{
    // defaults
    if (args()<4) d = missingof(to)
    
    // consistency checks
    if (eltype(x)!=eltype(from))  _error(3250, "{it:x} and {it:from} must be of same type")
    if (length(from)!=length(to)) _error(3200, "{it:from} and {it:to} must have same length")
    if (eltype(to)!=eltype(d))    _error(3250, "{it:to} and {it:d} must be of same type")
    if (length(d)!=1) {
        if (length(d)!=length(x)) _error(3200, "{it:x} and {it:d} not conformable")
    }
    
    // nothing to translate
    // - empty x
    if (!length(x)) return(J(rows(x), cols(x), missingof(to)))
    // - empty dictionary
    if (!length(to)) {
        if (length(d)!=1) return(J(1,1, rows(x)==rows(d) ? d : d'))
        return(J(rows(x), cols(x), d))
    }

    // run hash algorithm
    return(_mm_crosswalk_hash(x, from, to, d))
}

transmorphic vector _mm_crosswalk_hash(
    transmorphic vector x,
    transmorphic vector from,
    transmorphic vector to,
    transmorphic vector d)
{
    real scalar  i
    transmorphic A
    transmorphic scalar a
    real matrix notfound
    transmorphic vector res
    
    // build dictionary
    A = asarray_create(eltype(from))
    for (i=length(from); i; i--) asarray(A, from[i], to[i])
    
    // case 1: d is scalar
    if (length(d)==1) {
        res = J(rows(x), cols(x), missingof(to))
        asarray_notfound(A, d)
        for (i=length(x); i; i--) res[i] = asarray(A, x[i])
        return(res)
    }
    
    // case 2: d is non-scalar
    res = rows(x)==rows(d) ? d : d'
    notfound = J(0,0,.)
    asarray_notfound(A, notfound)
    for (i=length(x); i; i--) {
        a = asarray(A, x[i])
        if (a==notfound) continue // using asarray_contains() would be slower
        res[i] = a
    }
    return(res)
}

end
